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PREFACE 


This  report  is  one  of  a  series  of  studies  on  the  behavior  of 
materials  under  stress  conducted  under  Air  Force  Project  RAND.  The 
investigation  was  initiated  at  a  time  when  composite  materials  had 
begun  to  make  substantial  impact  on  new  airframe  design  under  Air 
Force  sponsorship. 

An  important  aspect  of  studies  of  the  application  of  composite 
materials  is  knowledge  of  their  failure  mechanisms.  Theories  to  date 
treat  only  part  of  the  problem — elastoplastic  behavior  of  the  composite 
material — and  have  not  been  able  to  treat  the  actual  failure  mode  of 
the  composite.  The  actual  failures  that  have  been  observed  seem  to 
occur  with  the  initiation  of  microscopic  cracks  at  the  boundary  of 
the  filament  and  the  matrix  and  their  subsequent  propagation  through¬ 
out  the  composite.  The  methodology  presented  in  this  report  is  a  first 
attempt  to  realistically  model  this  failure.  The  report  complements 
other  ongoing  research  and  development  now  being  sponsored  by  the 
Air  Force,  and  should  also  be  useful  to  the  composite  design  commun¬ 
ity. 

Dr.  Adams,  an  Associate  Professor  of  Mechanical  Engineering  at 
the  University  of  Wyoming,  is  a  consultant  to  The  Rand  Corporation. 
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SUMMARY 


The  behavior  of  materials  under  mechanical  stress  can  be  divided 
into  three  distinct  regimes:  (1)  linear  elastic  response  up  to  the 
elastic  limit,  (2)  inelastic  behavior  beyond  the  elastic  limit  and  up 
to  that  loading  at  which  first  failure  occurs  locally,  and  (3)  subse¬ 
quent  crack  propagation  and  total  composite  failure.  The  crack  ini¬ 
tiation  and  its  subsequent  propagation,  described  above  as  the  third 
regime,  is  the  subject  of  the  present  study. 

A  method  has  been  developed  for  predicting  the  strength  of  a 
unidirectional  composite  material  in  terms  of  its  micromechanical  re¬ 
sponse  to  an  applied  stress.  It  includes  elastoplastic  material  be¬ 
havior,  local  failure  that  initiates  a  crack,  and  propagation  of  the 
crack  to  cause  total  failure  of  the  composite. 

Although  the  basic  methodology  is  applicable  to  the  general  prob¬ 
lem,  a  specific  loading  condition — transverse  normal  loading — has  been 
selected  for  detailed  analysis.  It  is  often  the  transverse  properties 
that  limit  the  performance  of  the  composite  system.  Therefore,  this 
loading  condition  is  of  particular  interest  to  composite  materials 
technology  because  of  the  inherently  low  transverse  strength  of  most 
high-performance  composites. 

The  basic  principles  of  the  theory  of  plasticity  have  been  com¬ 
bined  with  a  finite  element  numerical  analysis  technique.  The  result 
is  a  rigorous  analysis  procedure  capable  of  accurately  modeling  the 
complex  boundary-value  problem  being  considered,  an  initial  step  to¬ 
ward  the  ultimate  goal  of  accurately  predicting  the  strength  of  a 
material.  A  complete  digital  computer  program  has  been  developed  as 
part  of  the  investigation,  permitting  the  ready  application  of  the 
analysis  to  practical  engineering  problems.  Because  the  primary  goal 
of  the  study  was  to  develop  a  method  of  analysis  and  to  write  an  as¬ 
sociated  computer  program,  only  limited  numerical  results  have  been 
obtained  to  date.  These  are  discussed  in  detail,  and  examples  dem¬ 
onstrating  the  type  of  information  which  can  be  obtained  are  given. 

Finally,  suggestions  for  a  number  of  possible  refinements  and 
extensions  of  the  present  analysis  are  outlined,  to  encourage  the 
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reader  to  use  and  to  further  develop  the  very  promising  analysis  method 
discussed  in  this  report.  Problem  areas  for  research  range  from  im¬ 
proving  the  accuracy  of  the  basic  finite  element  solution  technique  to 
developing  a  more  representative  model  of  the  propagating  crack. 
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A  =  6t2Kw 

O  I 

=  area  of  nth  triangular  finite  element 
a,  b  =  dimensions  of  the  finite  element  grid 

B  -  3To  +  (1  +  v)t 

D  =  A  +  ESj 
E  =  Young's  modulus 
f^j  =  element  node  point  forces 

[H]  =  plane  strain  elastoplastic  stiffness  matrix 
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[K]  =  element  stiffness  matrix 

[K]  =  array  stiffness  matrix 

2M  =  tangent  modulus  of  the  T  -  curve 
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n  =  total  number  of  node  points  in  array 
s_  =  deviatoric  stress  components 
AT  =  incremental  temperature  change 

Ca6  =  sa3  +  Vs335a6 

=  displacement  components 

Au,  Av  =  boundary  displacement  increments 

=  sum  of  forces  acting  at  a  node  point 
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a  =  coefficient  of  thermal  expansion 


engineering  shear  strain  components 


u  =  ultimate  value 


n  =  normal 
t  =  tangential 
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I .  INTRODUCTION 


The  behavior  of  materials  can  be  divided  into  three  distinct  regimes: 
(1)  linear  elastic  response  up  to  the  elastic  limit,  (2)  inelastic  be¬ 
havior  beyond  the  elastic  limit  and  up  to  that  loading  at  which  first 
failure  occurs  locally,  and  (3)  subsequent  crack  propagation  and  total 
composite  failure.  Behavior  within  the  third  regime  is  the  subject  of 
this  report. 

This  investigation  is  a  direct  extension  of  work  previously  reported 
in  Ref.  1,  which  had  as  its  main  topic  behavior  within  the  second  of  the 

regimes  given  above.  Included  in  Ref.  1  is  a  brief  review  of  the  his¬ 
torical  development  of  micromechanical  analyses  as  applied  to  the  impor¬ 

tant  problem  of  understanding  the  transverse  normal  loading  behavior  of 
a  unidirectional  composite. 

Transverse  normal  loading  is  of  special  interest  in  relation  to  uni¬ 
directional  composites  because  of  the  very  low  transverse  stiffness  and 
strength  characteristics  of  such  materials  relative  to  their  stiffness 
and  strength  in  the  direction  of  reinforcement  orientation.  Since  com¬ 
plex  laminated  composite  systems  are  constructed  of  individual  unidirec- 
tionally  reinforced  laminae  bonded  together,  this  highly  anisotropic 
behavior  exists  within  the  individual  lamina  even  when  the  laminated 
system  is  constructed  to  be  nearly  isotropic  with  respect  to  its  gross 
properties.  It  is  often  the  transverse  properties  that  limit  the  per¬ 
formance  of  the  composite  system.  But,  since  these  properties  are  very 
low,  there  is  a  high  potential  for  improvement. 

For  example,  the  transverse  tensile  strength  of  a  unidirectionally 

2 

reinforced,  epoxy  matrix  composite  is  typically  only  6000  to  8000  lb/in.  , 

whereas  the  longitudinal  tensile  strength  is  often  150,000  to  200,000 
2 

lb/in.  or  more.  The  tensile  strength  of  the  epoxy  matrix  itself  is 

2 

typically  12,000  to  15,000  lb/in.  .  Thus,  with  respect  to  transverse 
tensile  strength,  the  presence  of  the  filaments  actually  degrades  the 
matrix  material. 

A  detailed  theoretical  analysis  is  not  needed  to  establish  these 
characteristics  of  existing  materials;  they  are  regularly  demonstrated 
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experimentally  in  the  process  of  setting  material  property  design  allow¬ 
ables.  Rather,  the  purpose  of  this  micromechanical  analysis,  as  well  as 
those  previously  developed,  is  to  provide  a  better  understanding  of  the 
complex  interactions  between  reinforcing  filaments  and  the  surrounding 
matrix  material.  With  such  an  analysis,  it  is  possible  to  perform  a 
systematic  study  of  the  many  variables  that  influence  material  behavior. 
Such  variables  include  filament  shape,  filament  packing  geometry  and 
spacings,  filament  and  matrix  mechanical  properties,  interface  charac¬ 
teristics,  the  addition  of  filament  coatings,  and  the  presence  of  voids. 
Each  of  these  can  be  accurately  modeled  using  an  analysis,  and  their 
influence  on  the  gross  behavior  of  the  material  system  can  be  isolated. 
These  same  variables  are  usually  difficult,  if  not  impossible,  to  con¬ 
trol  accurately  in  an  experiment. 

This  analysis  of  crack  initiation  and  the  subsequent  propagation 
that  leads  to  total  failure  is  particularly  directed  toward  the  goal  of 
predicting  the  ultimate  strength  of  a  composite  material.  There  is 
currently  no  satisfactory  method  of  prediction.  However,  since  the 
strength  of  a  material  is  governed  by  local  effects  (as  opposed  to  com¬ 
posite  stiffness,  which  is  a  gross  response  averaged  over  the  entire 
material)  ,  it  is  reasonable  to  expect  that  a  local  or  micromechanical 
analysis  will  be  required.  This  investigation  represents  an  initial 
step  toward  the  ultimate  goal  of  accurately  predicting  the  strength  of 
a  material.  The  analysis  is  based  firmly  on  the  principles  of  solid 
mechanics.  However,  the  specific  methodology  for  modeling  a  crack  and 
its  propagation  is  relatively  crude.  In  subsequent  work,  methods  of 
applying  the  principles  of  fracture  mechanics,  which  are  rapidly  being 
developed  and  improved  upon,  should  be  introduced.  The  finite  element 
methodology  that  is  used  to  obtain  the  numerical  results  given  in  Sec.  IV 
also  is  being  developed  very  rapidly  at  the  present  time.  Improved 
element  representations  and  larger ,  more  economical  digital  computer 
facilities  will  have  a  very  favorable  influence  on  the  subsequent  pro¬ 
gress  of  micromechanical  analyses  of  this  type. 

A  brief  discussion  of  the  major  computational  aspects  of  the 
finite  element  computer  program  and  the  data— plotting  computer  program 
that  were  developed  as  part  of  this  investigation  appears  in  Sec.  III. 


To  the  author’s  knowledge,  this  investigation  represents  the  first 
attempt  to  predict  crack  initiation  and  propagation  in  a  composite 
material. 
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II.  METHOD  OF  ANALYSIS 


PREVIOUS  INVESTIGATIONS 

Reference  1  gives  details  of  the  author's  elastoplastic  analysis 
of  material  behavior  beyond  the  elastic  limit  and  up  to  that  trans¬ 
verse  normal  loading  at  which  first  failure  occurs  locally.  The  tangent 
modulus  method  of  treating  nonlinear  material  behavior  and  the  Prandtl- 
Reuss  flow  rule  are  employed.  A  plane  strain  condition  is  assumed,  and 
constant  strain  triangular  finite  elements  are  used.  The  stress-strain 
relationships  representing  each  of  the  constituent  materials  are  input 
point  by  point  so  that  any  monotonically  increasing  stress-strain  curve 
can  be  considered;  actual  experimental  data  are  typically  used.  Normal 
displacements  at  the  boundaries  must  be  specified  in  order  to  satisfy 
symmetry  conditions  for  the  boundary-value  problem  to  be  solved.  It  is 
usually  of  greater  practical  interest,  however,  to  be  able  to  specify 
instead  the  average  normal  stresses  on  each  boundary  (or  more  precisely, 
the  ratio  of  these  average  normal  stresses) .  To  maintain  this  desired 
ratio  exactly  for  each  solution  increment  requires  solving  two  separate 

boundary-value  problems  and  superimposing  the  results  in  the  proper  ratio. 

(2) 

This  has  been  a  standard  procedure  in  prior  works.  To  reduce  the  re¬ 

quired  computer  solution  time  by  approximately  one-half,  a  procedure  was 
developed  in  Ref.  1  to  estimate  the  boundary  normal  displacements  re¬ 
quired  to  produce  the  desired  increments  of  average  boundary  normal 
stresses.  Any  deviation  from  the  desired  stress  ratio  was  then  corrected 
for  in  the  succeeding  increment.  Quite  satisfactory  results  were  obtained 
using  this  predictor-corrector  approach.  A  separate  computer  program  was 
also  written  to  prepare  the  output  data  for  machine  plotting  of  contours. 

All  of  these  features,  with  the  exception  of  the  predictor-corrector 
method,  are  retained  in  the  analysis  of  crack  initiation  and  propagation 
reported  here.  The  predictor-corrector  method  was  found  to  be  insuf¬ 
ficiently  sensitive  for  the  current  application,  as  will  be  discussed  in 
more  detail  in  a  subsequent  subsection. 
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FILAMENT  PACKING  GEOMETRY 

The  physical  problem  to  be  considered  is  that  of  transverse  normal 
loading  of  a  unidirectionally  reinforced  composite,  as  illustrated  in 
Fig.  la.  Assuming  that  the  reinforcing  filaments  are  oriented  in  the 
direction  of  the  z-axis ,  as  indicated,  the  x-y  plane  then  represents 
the  transverse  loading  plane.  A  normal  loading  component  in  the  x- 
direction  is  illustrated  in  Fig.  la. 

In  most  composite  materials,  the  filament  packing  is  either  com¬ 
pletely  random  or  exhibits  only  a  very  slight  amount  of  preferred  orien¬ 
tation,  e.g.,  the  nesting  of  filaments  as  the  filament  volume  is  increased. 
One  exception  is  the  relatively  large  diameter  filaments  such  as  those  of 
boron  (which,  at  approximately  0.004  in.  in  diameter,  are  an  order  of 
magnitude  larger  than  most  glass  and  graphite  filaments  in  common  use). 
These  large  diameter  filaments  are  processed  as  monolayer  tapes  and  tend 
to  retain  this  layering  characteristic  in  the  finished  part.  This  is 
particularly  true  when  a  woven-glass  scrim-cloth  backing  layer  is  used 
to  improve  the  manageability  of  the  tape  in  fabrication. 

A  rigorous  continuum  mechanics  analysis  requires  the  establishment 
of  a  specific  filament  packing  array,  be  it  a  periodic  regular  array  or 
a  random  one.  The  problem  of  random  filament  packing,  and  a  method  of 
analyzing  its  influence  on  material  properties,  was  discussed  in  Refs.  3 
and  4.  The  necessity  (in  terms  of  computer  storage  capacity  and  running 
time)  of  selecting  a  representative  material  region  sufficiently  small 
to  permit  a  detailed  finite  element  representation  and  yet  be  physically 
meaningful  makes  the  accurate  modeling  of  a  random  array  somewhat  more 
difficult.  Thus,  in  the  current  investigation,  where  the  primary  emphasis 
is  on  developing  the  concept  of  crack  initiation  and  propagation,  only 
regular  packing  geometries  will  be  considered.  The  influence  of  random¬ 
ness  can  be  incorporated  in  subsequent  investigations,  as  desired. 

Two  commonly  assumed  periodic  regular  arrays  of  filaments  (cross 
sections  in  the  x-y  plane)  are  indicated  in  Figs,  lb  and  1c.  The  rec¬ 
tangular  array  in  Fig.  lb  includes  the  square  array  as  a  special  case 
(when  b  =  a) .  The  diamond  array  (sometimes  called  a  staggered  array) 
in  Fig.  lc  includes  the  hexagonal  array  as  a  special  case  (when  b  =  /3*  a). 


(a)  Transverse  normal  loading  of  a  unidirectional ly 
reinforced  composite 
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b)  Rectangular  array 


(c)  Diamond  array 


Fig.l  —  Filament  packing  geometries 
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Because  of  the  assumed  periodicity,  a  typical  repeating  unit  can 
be  isolated  in  each  case,  as  indicated  in  Figs,  lb  and  lc.  Only  one 
quadrant  of  this  repeating  unit  must  be  analyzed  (two  options  are  shown 
in  the  case  of  the  diamond  array) ,  because  of  the  assumed  symmetry  of 
the  filament  packing  and  the  geometry  of  the  individual  filaments  about 
both  the  x  and  y  axes.  The  boundary  conditions  on  the  sides  of  the 
quadrant  in  Fig.  lb  and  the  quadrant  on  the  left  side  of  Fig.  lc  are 
identical.  The  sloping  boundary  of  the  element  shown  on  the  right  side 
of  Fig.  lc  (solid  lines)  is  defined  as  passing  through  the  center  of 
that  quadrant  of  the  repeating  unit,  i.e.,  it  divides  the  quadrant  in 
half.  Either  of  the  elements  of  Fig.  lc  can  be  used,  yielding  identical 
results.  The  element  on  the  left  involves  simpler  boundary  conditions; 
however,  that  on  the  right  contains  only  one-half  as  much  area  and  only 
one  segment  of  filament-matrix  interface.  Thus,  only  one-half  as  many 
finite  elements  are  required  to  define  the  area  with  the  same  resolution. 
This  makes  the  use  of  the  element  on  the  right  of  Fig.  lc  almost  man¬ 
datory,  based  on  economics  alone.  But  the  boundary  conditions  on  the 
sloping  boundary  are  more  complicated,  particularly  for  an  elastoplas tic 
analysis.  This  will  be  discussed  in  greater  detail  in  Sec.  VI.  Only 
rectangular  arrays  will  be  analyzed  in  detail  here  since,  to  repeat,  the 
main  emphasis  is  on  developing  the  concept  of  crack  initiation  and 
propagation. 

The  first  quadrant  of  a  typical  repeating  element,  as  indicated  in 
Fig.  lb,  is  shown  in  Fig.  2.  The  cross-sectional  geometry  of  the  fila¬ 
ment  is  arbitrary  within  the  restriction  that  it  must  be  symmetrical 
about  both  the  x  and  y  axes.  Thus,  a  filament  region  of  any  arbitrary 
shape  drawn  in  the  first  quadrant  (Fig.  2)  is  allowable. 

At  one  time  there  was  considerable  interest  in  using  (glass)  fila¬ 
ments  of  special  shapes  to  enhance  specific  composite  properties.  For 
example,  rectangular  or  hexagonal  cross  sections  permit  very  high  packing 
densities  (theoretically,  filament  volumes  as  close  to  100  percent  as 
desired) .  Circular  filaments  in  square  and  hexagonal  arrays  are  limited 
to  maximum  filament  volumes  of  78.5  and  90.7  percent,  respectively  (when 
filaments  are  touching  each  other).  Thus,  the  composite  axial  strength 
and  stiffness  will  be  higher  for  the  specially  shaped  filaments.  However, 
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the  difficulty  of  accurately  placing  the  individual  small  filaments 
during  fabrication  to  achieve  the  high  packing  densities,  the  fact  that 
other  important  composite  properties  such  as  transverse  normal  and 
longitudinal  shear  strength  are  typically  degraded,  and  perhaps  most 
importantly,  the  introduction  of  high  modulus  filaments  such  as  boron 
and  graphite  that  are  not  produced  by  drawing  from  the  melt  (their 
shapes  being  relatively  fixed)  have  resulted  in  a  general  loss  of  inter¬ 
est  in  shaped  filaments  as  a  method  of  controlling  composite  properties. 

Boron  filaments  (boron  deposited  on  a  tungsten  substrate)  are  cir¬ 
cular  in  cross  section.  Graphite  filaments  produced  from  a  PAN  (poly¬ 
acrylonitrile)  precursor,  which  is  circular,  are  approximately  circular 
also.  On  the  other  hand,  graphite  filaments  produced  from  a  rayon  pre¬ 
cursor  retain  the  irregular  (crenulated)  cross-sectional  shape  of  this 
precursor  material.  Thus,  there  may  be  a  renewed  interest  in  the  future 
in  analyzing  the  influence  of  noncircular  filament  cross  sections.  The 
composite  longitudinal  shear  strengths  resulting  from  the  use  of  rayon 
precursor  graphite  filaments,  however,  are  typically  low  relative  to 
PAN  precursor  graphite  filament-reinforced  composites.  Whether  this  is 
due  to  the  irregular  shape  of  the  rayon  precursor  filaments  (creating 
unfavorable  local  shear  stress  concentrations)  or  to  some  other  effect 
has  not  been  established.  This  is  a  potential  future  application  for 
an  analysis  similar  to  the  one  reported  here.  In  the  meantime,  since 
this  problem  of  low  shear  strength  has  not  been  resolved,  more  and  more 
emphasis  is  being  placed  on  PAN  (and  other)  precursors,  all  of  which 
result  in  a  graphite  filament  of  approximately  circular  cross  section. 

An  interesting  new  organic  filament  developed  by  the  DuPont  Co.,^ 
designated  PRD-49-I,  has  very  attractive  dielectric  as  well  as  mechanical 
properties  and  a  low  density.  It  is  also  circular  in  cross  section. 

Thus,  current  interest  in  analyzing  the  influence  of  noncircular 
filaments  is  not  very  great.  In  this  investigation,  only  examples  of 
circular  filaments  have  been  used,  but  it  should  be  emphasized  that  the 
analysis  procedure  is  no  more  difficult  for  filaments  of  noncircular 
cross  section.  The  necessary  triangular  element  grid  can  be  constructed 
just  as  readily,  and  the  stress  and  displacement  continuity  requirements 
across  the  interface  are  automatically  satisfied  by  the  finite  element 
formulation,  independent  of  interface  shape. 
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BQUNDARY  CONDITIONS 

Specific  boundary  conditions  must  be  established  along  the  sides 

of  the  material  region  to  be  analyzed — the  first  quadrant  of  the  typical 

repeating  element  indicated  in  Fig,  2.  Since  a  two-dimensional  (plane 

strain)  analysis  is  being  used,  two  conditions  must  be  prescribed  at 

each  point  on  the  boundaries.  When  the  entire  composite  body  is  subjected 

to  any  combination  of  average  transverse  normal  stresses  a  and  a  ,  as 

x  y 

indicated  in  Fig.  1,  a  complex  state  of  stress  is  induced  locally  within 
the  material.  This  is  the  result  of  the  dissimilar  material  properties 
of  the  filaments  and  surrounding  matrix  and  also  a  result  of  interactions 
between  filaments.  Thus,  the  stress  distributions  along  the  boundaries 
of  the  first  quadrant  are  unknown  even  though  the  average  of  the  normal 
stresses  acting  along  each  side  must  equal  the  corresponding  applied 
stress  or  a  ,  from  equilibrium  considerations;  and  it  is  not  possible 
to  establish  normal  stress  boundary  conditions.  Because  of  the  assumed 
symmetry,  however,  the  rectangular  repeating  element  and  its  first  quad¬ 
rant  (bounded  by  lines  of  symmetry)  remain  rectangular  when  transverse 
normal  loads  are  applied;  the  normal  component  of  displacement  of  each 
point  on  a  given  boundary  of  the  first  quadrant  is  identical.  Thus,  an 
arbitrary  normal  displacement  of  each  boundary  can  be  specified,  and, 
because  of  the  assumed  symmetry,  no  shear  stresses  exist  along  the  boun¬ 
daries  of  the  first  quadrant. 

These  then  are  the  two  conditions  to  be  specified  at  each  boundary 
point:  a  uniform  normal  displacement  and  a  zero  shear  stress. 

Two  boundary  conditions  must  also  be  specified  at  each  point  on  the 
interface  between  dissimilar  materials,  e.g.,  at  the  filament-matrix 
interface.  Normally,  it  is  desired  to  completely  bond  the  filaments  to 
the  surrounding  matrix  during  the  fabrication  process.  In  the  analytical 
model,  this  corresponds  to  the  specification  of  interface  continuity  con¬ 
ditions  for  both  displacements  and  stresses,  i.e., 

f  m  f  m 

u  =  u  ,  u  =  u^ 


and 
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where  the  superscripts  f  and  m  represent  filament  and  matrix,  the  sub¬ 
scripts  n  and  t  represent  directions  normal  and  tangent  to  the  interface, 
and  a  and  T  represent  the  normal  stress  and  shear  stress.  This  is  usually 
termed  a  perfect  bond.  When  using  the  finite  element  method,  these  inter¬ 
face  continuity  conditions  are  automatically  satisfied  during  the  process 
of  assembling  elements  at  an  interface  node  point. 

In  certain  applications,  interface  boundary  conditions  other  than 
the  above  continuity  conditions  may  be  desired.  For  example,  it  may  be 
desired  to  model  the  effect  of  a  lack  of  interface  bond  in  a  local  region 
(perhaps  due  to  a  poor  fabrication  technique)  or  some  similar  initial 
defect.  The  interface  boundary  conditions  can  be  modified  accordingly. 

CONSTITUTIVE  RELATIONS 

The  material  region  of  interest,  as  introduced  in  the  preceding  sub¬ 
sections,  is  to  be  incrementally  loaded  along  its  boundaries.  As  will 
be  discussed  quantitatively  in  Sec.  IV,  a  suitable  matrix  material  for 
use  in  a  composite  is  typically  selected  for  its  ability  to  deform  ex¬ 
tensively  without  failing.  This  permits  the  large  strain  concentrations 
induced  in  the  matrix  region  between  adjacent  filaments  to  be  accommodated 
more  readily.  In  real  materials,  most  of  this  strain-to-f ailure  occurs 
beyond  the  linearly  elastic  range  of  behavior.  Thus,  for  loading  beyond 
the  elastic  limit,  the  elas toplastic  response  of  the  material  must  be 
considered  and  a  suitable  constitutive  relation  derived. 

The  present  analysis  is  not  restricted  to  monotonic  loading.  Thus, 
material  unloading  from  the  elastoplastic  region  may  occur,  requiring 
that  an  unloading  criterion  be  established  also.  It  will  be  assumed  that 
the  material  unloads  linearly  elastically,  with  a  modulus  equal  to  the 
initial  Young’s  modulus.  If  reloading  subsequently  occurs,  it  is  assumed 
that  the  material  follows  this  same  linear  stress-strain  response  until 
it  attains  the  state  from  which  it  began  to  unload.  For  additional  load¬ 
ing  beyond  this  point,  the  response  is  assumed  to  continue  along  the 
original  elastoplastic  curve. 

The  general  constitutive  relations  that  govern  the  elastoplastic 
response  can  be  derived  as  follows.  For  each  loading  increment  beyond 
the  elastic  limit,  the  corresponding  increment  of  strain  at  any  point  in 


-12- 


the  material  region  can  be  separated  into  an  elastic  (recoverable)  and 
a  plastic  (irrecoverable)  portion,  i.e., 


e.  . 

ij 


+  e.  . 


(1) 


The  elastic  portion  follows  the  generalized  Hooke’s  law, 


(e)  _  1-2V 


•  c  ,  1+V  • 

an5,,  +  ^  s  . 

kk  ij  E  ij 


where 


6.. 
kk  ij 


(3) 


is  the  deviatoric  component  of  the  rate  of  stress  tensor  ,  6  „  is 
the  Kronecker  delta,  and  E  and  V  are  the  material  Young’s  modulus  and 
Poisson’s  ratio,  respectively.  The  usual  conventions  of  indicial  nota¬ 
tion  apply.  Latin  indices  have  the  range  1,2,  3  and  Greek  1,  2, 

The  plastic  portion  is  assumed  to  obey  the  Prandtl-Reuss  flow  rule: 


e  . (p)  =  As . .  (4) 

il  ij 


where  X  is  a  positive  scalar  function,  in  general  a  function  of  both 

time  and  the  spatial  coordinates,  and  s^  bears  the  same  relationship 

to  O as  s  „  .  to  a..,  indicated  in  Eq .  (3).  The  key  assumption  is  that 
ij  ij  13 

the  rate  of  change  of  the  plastic  strain  is  at  any  instant  proportional 
to  the  deviatoric  stress.  It  is  also  assumed,  as  in  most  plasticity 
theories,  that  there  is  no  permanent  change  of  volume.  In  other  words, 
the  plastic  component  of  the  mean  normal  strain  is  zero;  the  deviatoric 
components  of  the  plastic  strain  are  identical  to  the  plastic  strain 
components.  The  Prandtl-Reuss  flow  rule  is  particularly  applicable  to 
the  case  of  contained  plastic  flow,  which  is  typical  of  the  condition 
existing  in  the  composite  material  geometries  being  considered  here. 
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In  such  a  case,  the  region  of  elastoplast ic  material  behavior  is  typi¬ 
cally  surrounded  or  contained  by  elastic  material.  Thus,  the  magnitudes 
of  the  elastic  and  plastic  components  of  strain  are  comparable. 

The  Prandtl-Reuss  flow  rule  does  not  model  viscosity,  i.e.,  time- 
dependent,  effects.  However,  most  of  the  composite  materials  in  common 
use,  e.g.,  the  metal  matrix  composites  and  those  incorporating  a  high 
polymer  matrix  such  as  an  epoxy,  do  not  exhibit  a  significant  time- 
dependent  response  at  ambient  temperatures.  Thus,  this  is  not  usually 
an  important  limitation. 

The  Prandtl-Reuss  flow  rule  is  the  one  most  frequently  used  today 
primarily  because  of  its  more  general  applicability  and  ease  of  implemen¬ 
tation.  Other  flow  rules  that  incorporate  specific  restrictive  assump¬ 
tions  also  are  occasionally  used.  For  example,  the  assumption  of  complete 
incompressibility  may  greatly  simplify  an  analysis,  and  provide  adequate 
accuracy  in  some  problems.  Also,  if  the  plastic  flow  is  unrestricted 
rather  than  contained,  so  that  the  plastic  strains  are  much  larger  than 
the  elastic  strains,  the  latter  can  often  be  neglected,  resulting  in  a 
simpler  rigid-plastic  analysis. 

In  the  applications  considered  in  this  report,  however,  these  various 
simplifying  assumptions  are  not  adequate  and  will  not  be  used.  Rather, 
the  Prandtl-Reuss  flow  rule  as  represented  by  Eq .  (4)  will  be  incorporated 
into  the  analysis  directly. 

To  put  X  into  a  usable  form  for  present  purposes,  Eq .  (4)  can  be 
multiplied  by  itself  and  solved  for  X  to  give 


a.  Me, 


JJ 


(p))i 


(sklskl)E 


(5) 


Defining  the  octahedral  plastic  shear  strain  rate  as 


e  <p)  -  (i 

o  3  ij  ij 


and  the  octahedral  shear  stress  as 


T 

O 


s . .s  .  .) 
1J 


2 


(6) 


(7) 
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Eq .  (5)  can  be  rewritten  as 


thus  defining  the  positive  scalar  function  A  in  Eq .  (4),  which  becomes 


•  (P) 

•  (p)  o 

e.  .  =  -  s  .  . 

il  To  1J 


Now  a  relationship  must  be  established  between  and  tq.  This 

can  be  done  by  applying  Eq .  (9)  to  a  problem  for  which  both  and 

T  are  known.  For  example,  since  we  will  be  interested  in  using  actual 
o 

constituent  material  experimental  data  in  generating  numerical  results, 
we  can  employ  the  results  of  a  simple  uniaxial  test,  for  which 


all  ^  a22  °33  °12  013  a23  ° 


(10a) 


•  (p)  a  a  ;  (p)  _  ;  (p)  =  _  I  e  <p) 

e  t  U,  e 22  -  too  2  ^11 


(p)  =  ;  <P>  =  e  <p)  =  0 


F  =  E 

b12  13 


(10b) 


•  (p) 

Substituting  Eqs .  (10b)  into  Eq .  (6)  and  solving  for  gives 

;  <p)  -  Jik  (p) 


Since  and  are  uni^orin  throughout  the  test  specimen,  they  are 

functions  of  t  only,  and  Eq.  (11)  can  be  rewritten  as 


,  ,  d£  (p)  de  (p)  di  Jl  i 

A  (P)  =  Vo 2 _  =  J2  — - - -  =  ■  -  — 

£11  Vl  dt  dTQ  dt  2Mt 


(12) 
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where ,  by  definition,  21^  =  dTQ/d£o^P^  is  the  tangent  modulus  of  the 
octahedral  shear  stress-octahedral  plastic  shear  strain  curve  for  the 
elastoplas tic  material,  and 


s  .  .s  .  . 


T  = 
o  3x 


is  obtained  by  differentiating  Eq.  (7).  Substituting  the  rate  of  plastic 
strain  values  defined  by  Eqs.  (10b)  and  (12)  into  Eq,  (6),  we  obtain 

;  (p)  _  lo  „,x 


Substituting  Eq.  (14)  into  Eq .  (9)  gives 


(p)  Vii 

"  2xo«r 


This  is  the  form  of  the  Prandtl-Reuss  flow  rule  desired* 

Substituting  Eqs.  (2)  and  (15)  into  Eq,  (1),  the  constitutive  re¬ 
lations  can  be  expressed  as 

1  _  l“2v  2.  o  ,  1+V  *  .  ToSii 

£  .  .  —  G,  ,  6  .  ,  4*  — — —  s  .  ,  +  - ( 16 ) 

ij  3E  kk  ij  E  xj  2t0mt  v  j 

Substituting  Eqs.  (3)  and  (13)  into  Eq .  (16)  and  rearranging  terms  gives 


r  -  1+v  A  v  i  r  ,  SijSkiakl 

-  T  “U  -  E  Vhj  +  -fT 2 — 

o  T 


In  obtaining  the  last  term,  use  has  been  made  of  the  fact  that  no  plastic 
work  is  done  by  the  hydrostatic  component  of  the  applied  stress, i.e., 


SklSkl  ~  Skl(0kl  3  amm6kl)  "  skiakl 
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since 


Si  .a  =0 

kl  mm 

Equation  (17)  is  the  general  flow  rule  relating  the  total  strain 
and  stress  increments. 

Ideally,  a  full  three-dimensional  analysis  of  the  boundary-value 
problem  outlined  in  the  previous  subsection  would  be  preferred.  This 
would  permit  the  consideration  of  stress  and  strain  gradients  in  the 
direction  of  filament  reinforcement  (the  3-direction  indicated  in 
Fig.  la)  as  well  as  in  the  transverse  (1-2)  plane.  For  example,  it 
would  be  possible  to  analyze  the  behavior  of  discontinuous  reinforcements 
such  as  whiskers  and  chopped  fibers  (as  in  the  transfer  of  stresses  at 
the  fiber  ends  into  adjacent  fibers)  and  particulate-reinforced  com¬ 
posites.  It  would  also  facilitate  the  modeling  of  voids,  local  debond¬ 
ing,  and  similar  manufacturing  defects.  The  analytical  formulation  of 
a  truly  three-dimensional  model  is  straightforward,  and  the  above  deri¬ 
vation  is  directly  applicable.  However,  even  for  a  two-dimensional 
analysis,  the  complexity  of  the  boundary-value  problem  makes  a  numerical 
solution  technique  mandatory.  Large  amounts  of  digital  computer  storage 
are  required  and  the  computer  running  times  are  relatively  long.  The 
present  generation  of  digital  computers,  e.g.,  the  IBM  Model  360  and  its 
counterparts,  do  not  have  sufficient  storage  capacity  and  computational 
speed  to  make  practical  three-dimensional  analyses  feasible.  The  same 
type  of  comments  applied  to  the  analysis  being  presented  here  just  a 
few  years  ago,  when  computers  such  as  the  IBM  Model  7094  and  its  counter¬ 
parts  were  typically  being  used — it  would  have  been  difficult  to  over¬ 
come  storage  limitations,  and  computer  running  times  required  would  have 
been  prohibitive  in  most  cases.  In  the  future,  solutions  of  truly  three- 
dimensional  practical  problems  will  undoubtedly  be  routinely  performed. 
Now,  however,  certain  idealizations  must  be  made  in  order  to  make  the 
problem  amenable  to  practical  solution. 

One  way  of  achieving  this  is  to  reduce  the  problem  to  a  two- 
dimensional  formulation.  Since  we  are  interested  in  composite  material 
response  in  the  transverse  plane,  it  is  logical  to  assume  that  stress 
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and  strain  gradients  in  the  axial  direction  (the  3-direction  in  Fig.  la) 
are  negligible  relative  to  those  in  the  transverse  plane.  This  is  the 
classical  plane  stress  or  plane  strain  condition.  A  plane  stress  con¬ 
dition  assumes  that  the  stress  in  the  direction  normal  to  the  plane  of 
interest  (the  component  of  stress,  normal  to  the  1-2  plane  in 
Fig.  la)  is  zero,  and  is  usually  associated  with  the  analysis  of  thin 
plates  (thin  in  the  3-direction)  subjected  to  in-plane  loadings.  But 
for  a  continuous-filament-reinforced  lamina,  as  indicated  in  Fig.  la, 
the  dimension  in  the  direction  of  reinforcement  (the  3-direction)  is 
typically  very  large.  This  corresponds  more  closely  to  the  plane  strain 
condition,  i.e.,  zero  displacement  in  the  3-direction.  Since  there  are 
assumed  to  be  no  gradients  in  the  3-direction,  the  strains  and  cor¬ 
respondingly  the  increments  of  strain  e^3  are  equal  to  zero,  i.e., 


9u -  9u« 

c  =  — i  4-  — ~  -  Q 

e13  9x3  9xx  U 


9u  9u« 

e  =  — L  4 — _  =  Q 
b23  9x3  9x2 


(18) 


9u, 

:33  =  dZ 


=  0 


where  u^  and  x^  represent  the  displacement  increment  components  and 
spatial  coordinates,  respectively. 

For  plane  strain, 


“l3  =  °23  =  °  (19) 

as  can  be  verified  by  substitution  into  Eq.  (17),  using  Eqs .  (18). 

Although,  by  definition,  the  axial  strain  increment  £^3  is  zero 
everywhere,  the  axial  stress  increment  0^3  is  not.  It  is  not  an  inde¬ 
pendent  quantity,  however,  and  can  be  expressed  in  terms  of  the  other 
nonzero  stress  increment  components  as  follows.  Setting  i  =  j  =  3  in 
Eq.  (17),  we  obtain 
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(1  +  V)  a33  -  vakk  +  A  s33skiakl  0 


where  A  =  6x^11^.  Solving  Eq .  (20)  for  ^335 


^33  Wyy  +  A  S33S33^33  +  A  s33sy6^y6  ° 


.  /  A  4*  Es33s33  .  e 

a33  - A -  *  WYY  '  A  s33sy«°y6 


where 


°33  "  D  (AWYY  “  ES33SY«V 


D  =  A  +  Es^ 

and  where,  as  previously  noted,  the  Greek  indices  have  the  range  1,  2. 
Utilizing  Eqs .  (19),  Eq .  (21)  can  be  substituted  into  Eq.  (17)  to 


, ,  /  Ava  -  Es00s  »a,  . 

•  =  1+v  •  _  V  -  •  YY  33  Y<S  Yd 

aB  E  aB  E  aB  YY  D 


s  ns  »a  r  s  „s„_  /  Ava, 
aB  Y^  Y<S  ,  aB  33  3 

A  A 


-  Es~-s...x.a_ 


Collecting  terms  in  a  and  gives 


1+V  •  _  V  .  ,D  +  AVs  VSaBS33  • 

eaB  E  aaB  E  °a8'  D  ’  D  YY 


'  v  *  ,  SaBSY|S  E  2  • 

+  D  (Sa6s33SY6  +  A  AD  saQs33sy6  °y 6 


or,  making  use  of  the  relation  D  =  A  +  Es.^, 


*  _  1+v  *  V  „  A(1  +  v)  V<^aBS33 

aB  E  aaB  ”  E  °a8  D  D 


xBS  33  ”1  • 

“D—  J  a. 


+  v  6  s  s  +  .S.aBsl6 
D  °aB  33  y«  / 


^  +  ^33  "  ^33 


Removing  the  common  denominator  D  and  rearranging  terms  gives 


^aB  "  "IT  V  +  h  'V(1  +  V)  I  6aBaYY 


+  vsQ0 (s  0  -  s 006  0)a 
33  a8  33  aB  YY 


+  (vs006  0  +  s  n)s  ra  r 

33  aB  aB  Y<$ 


Equation  (22)  is  equivalent  to  Eq .  (9)  of  Ref.  4,  as  demonstrated  by 
Swedlow  as  follows.  Let 


saB  taB  Vs336aB 


Equation  (22)  becomes 


®a6  "  TT  “as  +  5  -v<1  +  v)  I  'Vrr 


+  vs33  CaB  -  (1  +  V)s33SaB  a-f,  +  'crtsVYS 


Collecting  terms  on  v(l  +  v)  and  t^g  gives 


eaB  E  °a8  +  D  v(1  +  v)  (E  +  s33)<5aBaYY 


+  t  0(vs00a  +  s  m  r) 

aB  33  YY  y&  Y<5 
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Noting  that  (|  +  s^)  =  a  =  and»  as  previously  defined, 

(vs o  +  s  o)  =  t  o,  the  above  equation  reduces  to 

33  yo  yo  yo 


-v(l  +  v) 


^a6ayy  a6  y6  y6 


*  l+v  /*  *  ?  \  .  1  ,  .  * 

ea3  T  ^aaB  ^yy  aB  D  a8  y6  y6 


which  is  Eq.  (5)  of  Ref.  7. 

The  finite  element  numerical  analysis  will  be  formulated  and  solved 
in  terms  of  displacement  components.  Thus,  the  inverse  of  Eq.  (25)  is 
needed  so  that  stresses,  and  thus  node  point  forces,  can  be  readily  ex¬ 
pressed  in  terms  of  strain  components,  which  in  turn  are  simple  functions 


of  the  displacements. 

The  inverse  of  Eq .  (25)  is 


V  c 

ea8  l-2v  £YY°aB 


saBSY(S£'Yd 

3t2  +  (1  +  v) A/E 

o 


Methods  of  obtaining  this  inverse  are  given  in  Appendix  A.  (It  should 
be  noted  that  there  is  a  minor  typographical  error  in  the  inverse  given 
as  Eq.  (6)  in  Ref.  7;  MT  should  be  2^— the  factor  2  is  missing.) 

Equation  (26)  will  be  used  extensively  in  the  subsequent  development 
of  the  finite  element  representation.  For  this  purpose  it  is  convenient 
to  expand  Eq .  (26)  into  the  three  equations  it  represents,  i.e.,  CT^, 
a00,  and  a  10.  The  index  notation  will  also  be  changed  at  this  point 
to  conventional  engineering  notation.  Thus,  a^,  a22’  ai2’  £11*  e22’ 
e12  will  be  replaced  by  ax,  a  ,  Txy ,  £x,  £y,  \  Yxy,  respectively.  This 
notation  is  more  consistent  with  that  used  in  previously  published 
micromechanics  and  finite  element  analyses. 

Note  in  particular  that  £ ^  ^as  been  replaced  by  -j  Y  ,  t0  convert 
from  tensoral  shear  strain  to  engineering  shear  strain.  Engineering 
shear  strain  will  be  used  throughout  the  remainder  of  this  analysis. 
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Expanding  Eq .  (26): 


•  _  E  ,  v 

°x  1+v  (1  +  l-2v 


311X  •  ,  ,  v  slls22x  • 

T")ex  +  (l^  "  B  )£y 


2slls12.  1* 

(  B  >  2Yxy 


•  _  E  1— v  sll.  •  ,  ,  v  S11S22,. 

x  1+v  1-2V  “  B  )£x  +  l-2v  B  )£y 


,S11S12N • 

"  (_5— )YKy 


where 


B  -  3xJ  +  (1  +  «)  | 


Correspondingly,  for  and  t  , 


:  E  ,  V  slls22x •  .  ,  1— y  s22, • 

°y  1+V  (l-2v  B  )  x  +  (l-2v  B  )£y 
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Equations  (27)  ,  (28) ,  and  (29)  can  be  expressed  in  matrix  notation  as 


[a]  =  [H][e] 


(30) 


-22- 


where 


[a]  = 


—  - 

a 

€ 

X 

X 

. 

• 

* 

a 

,  [e]  = 

£ 

y 

y 

. 

♦ 

T 

xy 

xy 

[H]  = 


E 

1+V 


( 


1-v  11 


1-2V  B 


)  ( 


^l-2v  B 


S11S12 

( - — ) 

V  B 


V 


S11S22 


1-2V  B 

2 


S11S12 

)  (-  -J^) 


S  11s 22,  .  1-V  _  _2_2n 

)  \  1  O  „  v  T>  '  V 


l-2v  B 


s22s12, 

(  B  ' 


B 


S22S12, 


1  12 

(i  - 

k2  B  } 


Note  that  IH]  is  a  symmetric  matrix. 

FINITE  ELEMENT  REPRESENTATION 

Even  when  the  filament  packing  has  been  idealized  to  a  periodic 
regular  array  and  the  stress  distribution  reduced  to  a  two-dimensional 
plane  strain  condition,  the  resulting  elas toplas tic  boundary-value 
problem  is  not  amenable  to  a  closed-form  analytical  solution.  This  also 
was  true  in  prior  elastic  micromechanics  studies,  which  led  to  the  ap¬ 
plication  of  several  different  numerical  analyses  ,  including  stress 
function,  finite  difference,  and  finite  element  formulations.  An  his¬ 
torical  review  of  these  developments  may  be  found  in  Ref.  3.  Primarily 
because  of  the  greater  ease  of  representing  the  irregular  material  regions 
and  boundaries  (such  as  those  indicated  in  Figs.  1  and  2),  the  finite 
element  analysis  has  clearly  emerged  as  the  preferred  method. 

The  basic  finite  element  analysis  methodology  has  become  well- 
developed  in  recent  years,  and  has  been  used  extensively  in  structural 
analysis  applications.  A  number  of  texts  are  now  available,  e.g.,  that 
by  Zienkiewicz, ^  which  describe  the  method  in  detail.  Particular 
emphasis  today  is  on  the  development  of  improved  element  representations, 
including  three-dimensional  elements. 
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Constant  strain  triangular  elements  have  been  used  in  this  inves¬ 
tigation,  That  is,  the  displacement  field  in  each  triangular  element 
is  assumed  to  be  linear.  The  first  step  is  to  divide  the  material 
region  of  interest  into  triangular  elements,  as  is  indicated  in  Fig .  3. 

As  a  general  rule,  the  individual  elements  should  be  compact,  i.e., 
shaped  as  close  to  equilateral  triangles  as  is  practical  to  best  repre¬ 
sent  the  assumed  constant  strain  condition  within  each  element. 

The  amount  of  computer  solution  time  required  is  a  direct  function 
of  the  number  of  element  node  points.  The  strain  field  can  thus  be  most 
accurately  determined,  for  a  given  number  of  node  points,  by  using  small 
elements  in  regions  where  gradients  are  expected  to  be  high  (for  example, 
in  regions  between  directly  adjacent  filaments  such  as  along  the  x  and  y 
axes  of  Fig.  3),  and  larger  elements  elsewhere.  However,  this  procedure 
is  not  completely  satisfactory  for  the  present  crack  propagation  analysis 
As  the  numerical  results  will  indicate,  it  is  not  always  readily  apparent 
beforehand  where  the  crack  will  initiate,  or  where  it  will  propagate  to. 
Also,  once  the  crack  propagation  process  begins,  there  will  be  a  sig¬ 
nificant  redistribution  of  strains,  particularly  near  the  crack  tip.  To 
achieve  total  failure,  the  crack  must  eventually  propagate  across  the 
entire  material  region.  Thus,  it  is  better  to  keep  the  grid  array  rela¬ 
tively  uniform  in  size. 

Of  course,  one  constituent  material  may  be  much  stronger  than  the 
other  (or  others,  if  more  than  two  constituents  are  represented),  so 
that  crack  propagation  is  certain  to  occur  only  in  the  weaker  constituent 
In  this  case,  a  coarser  grid  can  be  used  in  the  region  of  strong  material 
in  order  to  conserve  node  points.  Such  a  grid  is  indicated  in  Fig.  3, 
where  the  circular  filament  region  is  represented  by  a  coarser  grid  than 
the  surrounding  matrix  region.  The  displacement  of  a  given  node  point 
is,  however,  a  function  of  the  displacements  of  all  of  the  immediately 
surrounding  node  points  (and  the  material  properties  of  the  elements 
they  bound).  Thus,  for  the  local  strain  gradients  to  be  reasonably 
represented  from  element  to  element,  it  is  desirable  to  use  a  transition 
from  small  to  large  elements,  as  is  indicated  in  Fig.  3,  rather  than  an 
abrupt  change.  Each  node  point  and  each  element  must  be  numbered  for 
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Fig.  3 —  Typical  finite  element  grid:  square  array,  40-percent 
filament  volume,  176  nodes,  304  elements 
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identification.  For  clarity,  in  Fig.  3  only  the  numbers  of  those  ele¬ 
ments  that  are  specifically  referred  to  in  Sec.  IV  or  Appendix  B  are 
shown. 

Having  constructed  a  suitable  finite  element  grid,  each  triangular 
element  is  defined  by  the  x,  y  coordinates  of  its  three  node  points, 
i,  j,  k,  as  indicated  in  Fig.  4 — i.e.  ,  x^ ,  y^ ,  x^  ,  y  ^  ,  x^,  y^.  Because 
of  the  nonlinear  (elastoplas tic)  material  behavior,  the  loading  is  to 
be  applied  in  small  increments  beyond  the  elastic  limit,  as  will  be 
described  in  the  next  subsection.  For  each  load  increment,  the  incre¬ 
ment  of  displacement  of  any  point  within  an  element  can  be  represented 
as 


u  =  a^x  +  a2y  + 
v  =  BjX  +  62Y  +  B3 


(31) 


where  u  and  v  are  the  x  and  y  components,  respectively,  of  the  incre¬ 
ment  of  displacement  of  any  point  in  the  element.  The  six  constants 
represented  by  the  a  and  8  coefficients  can  be  evaluated  by  substituting 
the  values  of  u,  v,  x,  and  y  for  each  of  the  three  node  points  of  the 
element  into  Eqs .  (31).  The  resulting  three  equations  containing  the 
three  unknown  a  constants  can  then  be  solved  simultaneously  to  evaluate 
,  a2 ,  Oi^.  The  remaining  three  equations  correspondingly  yield  8^  , 

B2 ,  8^.  Substituting  these  values  of  a  and  8  into  Eqs.  (31)  gives 

(Vk  -  Vj> +  yjkx  + 

+  [<vi  -  Vk>  +  >'kix  +  vj  "j 
+  [<xiyj  -  xjyi)  +  yijx  +  xjiy]  “kj 


(32) 


Fig.  4 — nth  triangular  element 


where 


v 


|[(xjyk  "  Vj}  +  yjkx  +  x 

(\yi  -  xiyk}  +  ykix  +  xiky 

(xiyj +  xjyi}  +  y±jx  +  xjiy 


A 

n 


■  i 


Vij> 


represents  the  area  of  the  triangular  element  n  and  the  notations 
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Xi  ~  Xj ’ 


y .  .  =  y .  -  y . 
ij  J3 


have  been  used.  These  displacement  increment  functions,  Eqs .  (32)  and 
(33),  assure  continuity  of  displacements  between  adjacent  elements, 
because  displacements  are  necessarily  equal  at  common  node  points  at 
each  end  of  the  common  side  and  vary  linearly  along  the  side. 

The  strain  increment-displacement  increment  relations 


3u  _9v 
3y  3x 


can  be  expressed  in  terms  of  the  node  point  displacement  increments  by 

taking  the  partial  derivatives  of  Eqs.  (32)  and  (33).  (Note  that,  as 

previously  discussed,  y  in  Eqs.  (34)  represents  the  increment  of  engi- 

xy 

neering  shear  strain,  and  not  the  tensoral  value.)  Equations  (34) 
become 


[E]  =  [6]  [6] 


where 
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Element  stress  increments  can  be  expressed  in  terms  of  node  point 
displacement  increments  by  substituting  Eq .  (35)  into  Eq.  (30): 

[O]  =  [H]  [9]  [5]  (36) 


These  stress  increments,  uniformly  distributed  in  the  element  and 
thus  along  its  boundaries,  can  be  replaced  by  statically  equivalent 
force  increments  acting  at  each  node  point  of  the  element.  These  equiv¬ 
alent  force  increments  will  be  defined  in  terms  of  their  x  and  y  co¬ 
ordinates  as 


[F] 


The  finite  element  analysis 
librium  of  the  forces  contributed  by  each  surrounding  element  at  each 
node  point  in  the  element  array.  Thus,  it  is  necessary  to  express  the 
node  point  force  increments  for  each  element  in  terms  of  the  node  point 
displacement  increments.  One  simple  procedure  is  to  impose  arbitrary 
(virtual)  node  point  displacements  and  to  equate  the  external  and  inter¬ 
nal  work  by  the  various  force  and  stress  increments  acting  on  the 


X. 

l 

Y. 
i 


X. 

J 


(37) 


is  based  on  the  establishment  of  equi- 


element . 
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In  general ,  such  forces  may  include  body  force  components  dis¬ 
tributed  over  the  element,  stress  distributions  due  to  temperature 
changes,  processing  residual  stresses,  and  so  forth.  These  effects 
can  be  introduced  in  a  straightforward  manner,  if  desired  (see,  for 
example,  Ref.  8  for  a  detailed  discussion) ;  they  are  not  included 
here . 

The  work  done  by  the  node  point  force  increments  (external  work) 
is  equal  to  the  sum  of  the  dot  products  of  the  individual  force  com¬ 
ponents  and  the  corresponding  displacement  increment  components,  i.e., 

WE  =  [6*]T  •  [F]  (38) 

where 


represents  the  virtual  displacements  at  the  element  node  points,  [F] 
is  defined  in  Eq .  (37),  and  the  superscript  T  denotes  the  transpose 
of  the  matrix. 

The  internal  work  per  unit  volume  done  by  the  stress  increments 
is 

W.  =  [e*]T  [a]  (39) 

where 

[£*]  =  [9]  [6*] 
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as  indicated  in  Eq .  (35);  [6*]  is  defined  in  Eq .  (38).  The  transpose 

•* 

of  [e  ]  is 

[e*]T  =  ([6]  [<S*])T  =  [5*]T  [0]T  (40) 

Substituting  Eqs .  (36)  and  (40)  into  Eq.  (39)  gives 

W.  =  [6*] T  [0]T  [H]  [9]  [6] 

The  total  internal  work  is  then  obtained  by  integrating  this  work  per 
unit  volume  over  the  volume  of  the  finite  element 

WI  =  /  [0]T  IH1  f0]  [h  d(vol)  (41) 

Equating  the  external  and  internal  work  done,  i.e.,  Eqs.  (38)  and  (41), 

[F]  =  [K]  [6]  (42) 

where 

[K]  =  f [0]T  [H]  [0]  d(vol) 

is  the  element  stiffness  matrix.  For  the  case  considered  here,  where 
the  stresses  and  strains  are  constant  throughout  the  element,  both  [0] 
and  [H]  are  constant  and  the  integration  is  trivial,  i.e., 

[K]  =  Ant0]T  [H]  [9]  (43) 

The  above  procedure  for  obtaining  the  required  relationship  between 
node  point  force  increments  and  node  point  displacement  increments, 

Eq.  (42),  is  applicable  for  any  assumed  strain  distribution  in  the  ele¬ 
ment.  However,  for  the  case  where  the  strains  (and  hence  stresses)  are 
assumed  to  be  constant  within  each  element,  an  even  simpler  and  perhaps 
more  physically  intuitive  procedure  can  be  followed,  as  was  done,  for 
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(9) 

example,  by  Foye.  The  triangular  element  of  Fig*  4  is  assumed  to 

be  inscribed  in  a  rectangular  region  of  material,  the  sides  of  which 
are  parallel  to  the  x-y  coordinate  axes.  The  stresses  throughout  this 
rectangular  region,  including  along  its  sides,  are  assumed  to  be  con¬ 
stant,  as  they  are  in  the  triangular  element.  This  rectangular  region 
can  then  be  divided  into  four  free  bodies,  as  indicated  in  Fig.  5;  the 
equivalent  force  increments  (per  unit  thickness  in  the  z-direction) 

acting  on  the  sides  are  as  shown  in  the  figure.  Since  the  stresses 

•  «  • 

O  ,  a  ,  T  are  constant  along  each  side,  the  equivalent  force  incre- 

^  y  .  ^y  . 

ments,  a  y..,  f..,  etc.,  can  be  assumed  to  be  acting  at  the  midpoints 
x  j  i  1 2 

of  the  sides,  as  indicated.  Summing  forces  in  the  x  and  y  directions 
for  each  of  the  three  free  bodies  surrounding  the  actual  finite  element 
gives 


x 


f .  .  =  a  y .  .  +  T  x., 
ij  xJ2 1  xy  ij 


*  y 

f .  .  =  a  x .  .  +  T  y  *  , 
ij  Y  iJ  xyJji 


x 


fjk  ■  Vkj  +  TxyXjk 

.  y  .  • 

f.:  -ax..  +  t  y,  , 
jk  y  jk  xy  kj 


(44) 


fik  =  Vik  +  TxyXki 


.  y  •  • 

f,:  =  a  x.  .  +  r  yM 
lk  y  Tci  xyyik 


These  are  the  statically  equivalent  force  increment  components  acting 
on  the  midpoints  of  the  sides  of  the  nth  triangular  element  of  area 
as  indicated  in  Fig.  5.  These  force  increments  can  then  be  divided 
equally  between  the  two  end  points  of  the  side  on  which  each  acts.  For 
example,  the  sum  of  the  x-components  of  these  divided  forces  at  the  node 
point  i  of  the  nth  element  is,  using  the  notation  of  Eq .  (37), 
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It  can  be  readily  shown  that  these  equations  are  equivalent  to 
Eq.  (42).  Integrating  Eq .  (39)  over  the  volume  of  the  nth  triangular 
element  gives 


W 


I 


[a]  d(vol) 


•  •k  T 

=  An[e  ]x  [a] 


(46) 


since  [s  ]  and  [a]  are  constant  throughout  the  element.  Substituting 

•  *  x 

Eq.  (40)  for  [e  ]  into  Eq .  (46)  and  equating  the  result  to  Eq .  (38) 
gives 


[F]  =  An[0]T  [a]  (47) 

where  [F]  ,  [0],  and  [a]  are  defined  in  Eqs .  (37),  (35),  and  (30),  re¬ 
spectively.  Equation  (45)  is  the  first  of  the  set  of  equations  repre¬ 
sented  by  Eq.  (47).  The  remaining  equations  are  similarly  equivalent. 
Equation  (36)  can  be  substituted  into  the  set  of  equations  represented 
by  Eq .  (45)  to  obtain  Eq .  (42). 

Equation  (47)  will  be  used  later  when  the  problem  of  accounting 
for  failed  elements  is  discussed. 

Equation  (42)  relates  the  node  point  force  increment  equivalents 
of  the  stress  increments  in  each  element  to  the  node  point  displacement 
increments.  If  body  forces  are  present  (a  situation  not  included  in 
this  investigation) ,  an  additional  set  of  node  point  force  increment 
equivalents  would  exist.  Also,  if  the  boundaries  of  the  material  region 
represented  by  the  finite  element  array  are  subjected  to  distributed 
external  loadings,  the  equivalent  force  increment  components  acting  on 
the  individual  boundary  nodes  must  be  included.  This  will  result  in 
an  additional  set  of  node  point  force  increments  [R] . 

Having  defined  all  of  the  equivalent  node  point  forces  acting  on 
each  element,  the  next  step  is  to  reassemble  the  individual  elements  in 
the  array.  The  net  force  at  each  node  point  in  the  array  is  then  ob¬ 
tained  by  summing  up  the  force  increment  components  contributed  by  the 
surrounding  elements : 
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[F]  =  [K]  [6] 


(48) 


where  the  bars  represent  summations  in  each  coordinate  direction  at 
every  node  point  in  the  array.  Thus,  [F]  and  [6]  will  be  column  vectors, 
each  with  2N  entries  ,  where  N  is  the  total  number  of  node  points  in  the 
array.  The  entries  in  [I?]  are  the  sum  of  the  components  of  force  at  each 
node  point  in  the  array  contributed  by  the  stresses  in  the  elements  sur¬ 
rounding  that  node  point,  i.e.,  the  sum  of  the  equivalent  force  components 
represented  by  Eq.  (42)  for  each  surrounding  element.  Because  each  node 
point  must  remain  in  equilibrium,  this  sum  must  equal  the  externally 
applied  force  component  at  that  node  point,  or  the  corresponding  entry 
in  [R] .  Thus,  for  interior  node  points  (where  no  boundary  force  com¬ 
ponents  are  acting,  and  the  entry  in  [R](is  therefore  zero),  the  sum  of 
these  equivalent  force  components  must  be  zero,  resulting  in  a  zero  entry 
in  [F] .  For  node  points  on  the  boundary,  the  sum  of  the  equivalent  force 
components  from  Eq.  (42)  must  equal  the  component  of  the  equivalent  boun¬ 
dary  force  acting  there,  the  corresponding  entry  in  [R] .  At  boundary 
node  points  where  an  equivalent  boundary  force  increment  component  of 
known  magnitude  is  applied,  the  sum  of  the  equivalent  force  increment 
components  due  to  element  stresses  must  equal  this  magnitude  in  order  to 
maintain  equilibrium,  i.e.,  the  corresponding  entry  in  [F]  must  equal 
this  known  value.  At  boundary  node  points  where  a  displacement  increment 
component  of  known  magnitude  is  applied,  the  associated  boundary  force 
component  increment  is  an  unknown  quantity,  to  be  evaluated  as  part  of 
the  solution  process.  The  corresponding  entry  in  [F]  is,  then,  an  un¬ 
known.  Thus,  the  2N  entries  in  the  column  vector  [F]  represent  the 
sum  of  the  x  and  y  components  of  the  equivalent  force  increments  at  the 
N  nodes  of  the  array  contributed  by  each  of  the  surrounding  elements. 

Since  most  of  the  nodes  in  a  typical  array  are  interior  nodes,  most  of 
the  entries  in  [F]  will  be  zero. 

The  [K]  in  Eq .  (48)  is  a  2N  x  2N  symmetric  positive  definite  matrix. 
Each  entry  in  [K]  is  the  sum  of  terms  arising  from  the  application  of 
Eq.  (42)  for  each  element  surrounding  a  given  node.  Since  the  total  num¬ 
ber  of  nodes  N  is  typically  large  relative  to  the  number  of  elements  (and 
therefore  the  number  of  nodes)  surrounding  a  given  node,  [K]  will  be  a 
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sparse  matrix,  i.e.,  most  of  its  entries  will  be  zero.  Advantage  will 
be  taken  of  this  sparseness  later  when  discussing  details  of  the  pro¬ 
cedure  for  solving  Eq .  (48)  for  the  unknown  displacement  increments. 

As  will  be  discussed  in  the  next  subsection,  the  boundary-value 
problems  to  be  solved  here  are  of  the  mixed-mixed  type;  at  most  of  the 
boundary  node  points,  a  force  increment  in  one  coordinate  direction  and 
a  displacement  increment  in  the  other  is  prescribed.  (Note  that  corre¬ 
sponding  components  of  force  and  displacement  cannot  both  be  prescribed.) 
Equation  (48)  can  be  rearranged  so  that  it  may  be  partitioned  in  the 
form 

(49) 


•  k  *  u 

where  [F_  ]  and  [_F  ]  are  the  known  and  unknown  node  point  force  incre- 

*  u  •  k 

ments ,  respectively,  and  [S_  ]  and  [j$  ]  are  the  corresponding  unknown 
and  known  node  point  displacement  increments.  The  unknown  force  incre- 
ments  [JF  ]  can  be  evaluated  directly  from  Eq .  (49)  once  the  unknown 
displacement  increments  [j$  ]  are 

[FU]  =  [K21] 

•  u 

[j5  ]  remains  to  be  solved. 

[Zk]  =  [Kn] 

which  can  be  rewritten  as 

[6U]  =  lKn]_1 


solved  for: 

tlU]  +  [K22]  [6k]  (50) 

From  Eq .  (49)  , 

[6U]  +  [K12]  [|k] 


[Fk3  -  [K12]  [6k]  (51) 
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Note  that  if  two  rows  of  [K]  are  interchanged,  the  corresponding 
columns  are  also  interchanged.  Thus,  the  symmetry  of  the  matrix  is 
preserved . 
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where  [K^]*"1  is  the  inverse  of  [K^]  .  In  practice,  the  challenging 
problem  is  to  develop  an  efficient  computational  procedure  for  obtain¬ 
ing  this  inverse,  since  the  order  of  [K.^]  typically  large* 

The  solution  of  Eq .  (51)  for  the  unknown  node  point  displacement 
increment  components  is  the  key  step  in  the  solution  procedure.  The 
unknown  force  increment  components  (at  node  points  where  the  corre¬ 
sponding  displacement  increment  components  are  known  initially)  can 
then  be  computed,  if  required,  from  Eq .  (50).  Strain  increment  com¬ 
ponents  within  each  finite  element  are  computed  by  substituting  the 
node  point  displacement  increments  into  Eq .  (35).  The  corresponding 
element  stress  increment  components  can  be  computed  either  by  substi¬ 
tuting  the  node  point  displacement  increment  components  into  Eq.  (36) 
or  by  substituting  the  computed  strain  increment  components  into  Eq. 

(30).  Thus,  the  complete  solution  is  obtained  for  the  given  applied 
stress  increment  at  the  boundary. 

The  detailed  computational  procedure  for  assembling  these  incre¬ 
mental  solutions  into  a  complete  analysis  can  now  be  introduced. 

INCREMENTAL  ANALYSIS  METHODOLOGY 

The  constitutive  relations  for  elastoplastic  material  behavior 
(which  were  derived  in  a  previous  subsection)  are  piecewise  linear  if 
the  slope  of  the  octahedral  shear  stress-octahedral  plastic  shear  strain 
curve  (the  tangent  modulus  2M^  in  Eq .  (12))  can  be  approximated  to  be 
constant  during  each  (small)  increment  of  loading.  Thus,  the  full  load¬ 
ing  of  a  material  well  into  the  nonlinear  (elastoplastic)  region  can  be 
treated  as  a  sequence  of  linear  loading  increments. 

Two  general  approaches  have  been  developed  for  treating  this  incre¬ 
mental  loading  problem:  the  method  of  initial  strains  and  the  tangent 
modulus  method.  Each  has  advantages  and  disadvantages  as  applied  to 
the  composite  materials  problem. 

The  method  of  initial  strains  is  based  on  the  fact  that  plastic 
strain  components  do  not  cause  a  change  in  stress.  At  the  beginning  of 
each  load  increment,  an  initial  guess  of  the  increments  of  plastic  strain 
that  will  occur  in  each  element  is  made.  (For  simplicity,  zero  values 
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may  be  assumed.)  Treating  these  values  as  initial  strains,  a  linear 
elastic  analysis  is  performed  for  the  prescribed  increment  of  boundary 
loading.  After  obtaining  the  resulting  stress  increments  for  each 
element,  a  revised  estimate  of  the  plastic  strain  increments  in  each 
element  can  be  computed  using  these  stress  increments,  the  material  con¬ 
stitutive  relations  (Eq.  (15)),  and  the  value  of  2MT  appropriate  to  the 
current  level  of  accumulated  stress  in  each  element.  These  revised 
estimates  of  the  plastic  strain  increments  are  then  used  as  initial 
strains  and  the  iterative  procedure  is  repeated  as  many  times  as  is 
necessary  to  achieve  adequate  convergence.  A  common  criterion  of  ade¬ 
quacy  is  that  the  current  estimates  differ  by  not  more  than  a  prescribed 
amount  from  those  obtained  in  the  previous  iteration. 

The  principal  advantage  of  the  method  of  initial  strains  is  asso¬ 
ciated  with  the  fact  that  only  linear  elastic  boundary-value  problems 
need  to  be  solved.  The  elastic  stiffness  matrix,  analogous  to  the 
elastoplastic  stiffness  matrix  [K]  in  Eq .  (48),  involves  only  the  elastic 
material  constants  E  and  V,  which  are  known  constants.  Thus,  the  stiff¬ 
ness  matrix  does  not  change  throughout  the  loading  process ,  which  means 
that  its  inverse  need  only  be  computed  once.  The  displacement  increments 
for  each  iteration  are  computed  directly  from  a  set  of  equations  anal¬ 
ogous  to  those  of  Eq.  (51).  Terms  similar  to  those  in  parentheses  in 
Eq.  (51)  change  from  iteration  to  iteration  and  increment  to  increment, 
but  the  stiffness  matrix  does  not. 

A  potential  disadvantage  of  the  method  of  initial  strains  is  the 
possible  slow  convergence  of  the  iteration  process.  In  practical  appli¬ 
cations  this  is  particularly  likely  to  occur  for  loadings  well  beyond 
the  elastic  limit,  where  the  slope  of  the  material  stress-strain  curve 
(the  tangent  modulus  2M^)  usually  becomes  small.  In  such  cases,  a  small 
increase  in  stress  in  the  element  corresponds  to  a  large  increment  of 
plastic  strain,  causing  the  iteration  process  to  converge  very  slowly, 
or  possibly  not  at  all.  Since  in  the  current  investigation  we  are 
interested  in  accurately  representing  the  material  response  within  each 
finite  element  to  full  failure  (crack  initiation  and  propagation) ,  this 
can  become  a  very  serious  problem. 


-38- 


The  tangent  modulus  method  of  analysis  of  elastoplastic  material 
behavior  uses  the  same  general  solution  techniques  that  were  developed 
previously  for  linear  elastic  analyses.  The  elastoplastic  constitutive 
relations,  Eq .  (17),  replace  the  simpler  linear  elastic  relations  (and 
reduce  to  them  as  a  special  case  when  the  plastic  strain  is  zero,  i.e., 
when  the  last  term  in  Eq .  (17)  is  zero,  corresponding  to  2Hj  =  °°)  .  This 
slope  of  the  governing  octahedral  shear  stress-octahedral  plastic  shear 
strain  curve  (21^  in  Eq .  (17)),  is  assumed  to  remain  constant  during 
each  load  increment .  A  new  value  is  determined  for  each  element  at  the 
beginning  of  each  load  increment.  Thus,  the  strain  increments  given 
by  Eq.  (17)  are  a  linear  function  of  the  stress  increments  during  each 
load  increment. 

The  principal  advantage  of  the  tangent  modulus  method  is  that  the 
solution  for  each  load  increment  is  obtained  directly  rather  than  by 
iteration.  However,  the  stiffness  matrix  [K]  in  Eq .  (48)  contains  a 
nonzero  value  of  2M^  for  each  element  that  has  exceeded  its  material 
elastic  limit,  and  these  values  usually  change  from  one  load  increment 
to  the  next.  Thus,  it  is  necessary  to  construct  a  new  [K]  at  the  begin¬ 
ning  of  each  load  increment,  and  invert  the  corresponding  submatrix  [K^] 
for  use  in  Eq.  (51).  It  is  particularly  the  inversion  of  [K^]  that  con¬ 
sumes  computer  time. 

In  actual  practice,  the  determination  of  which  method  will  require 
less  computation  time  will  depend  upon  the  rate  of  convergence  of  the 
method  of  initial  strains  solution.  For  example,  Marcal^10^  has  cited 
a  typical  problem  for  which  the  computing  time  required  to  achieve  equal 
accuracy  using  the  method  of  initial  strains  was  twice  as  long  as  for 
the  tangent  modulus  method.  A  similar  trend  can  be  expected  for  the 
composite  material  problem,  discussed  here,  since  large  strains  (to 
failure)  are  of  principal  interest. 

As  is  generally  true  when  employing  iterative  schemes,  it  is  very 
difficult  to  estimate  beforehand  the  amount  of  computation  time  required 
when  using  the  method  of  initial  strains.  Conversely ,  the  time  required 
per  load  increment  when  using  the  tangent  modulus  method  can  be  accurately 
estimated  in  advance,  because  it  is  directly  proportional  to  the  number 
of  unknown  node  point  displacements  to  be  solved  for  (and  the  bandwidth 
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of  the  diagonalized  stiffness  matrix  to  be  inverted,  as  will  be  discussed 
later)  . 

Considering  the  various  advantages  and  disadvantages  of  the  two 
approaches  discussed  in  this  subsection,  the  tangent  modulus  method,  as 
developed  for  use  in  the  previously  reported  inelastic  analysis, ^  has 
been  retained  for  this  investigation  of  crack  initiation  and  propagation. 
An  example  of  the  application  of  the  method  of  initial  strains  to  in¬ 
elastic  composite  material  behavior  has  been  given  by  Foye  and  Baker. 


INCREMENTAL  LOADING  TO  FIRST  FAILURE 

Elastoplas tic  material  behavior  beyond  the  elastic  limit  and  up  to 
that  transverse  normal  loading  at  which  the  stress  in  any  element  reaches 
its  ultimate  strength  value  is  the  subject  of  a  previous  report.  ^  Be¬ 
cause  the  problem  is  not  formulated  there  in  detail  and  because  a  number 
of  refinements  have  been  made  during  the  present  investigation,  a  fuller 
discussion  will  be  included  here. 

The  general  boundary  conditions  on  the  sides  of  the  rectangular 
first  quadrant  (the  region  over  which  the  solution  is  to  be  obtained; 
see  Fig.  2)  and  the  continuity  conditions  at  the  filament-matrix  inter¬ 
face  were  defined  in  a  previous  subsection.  A  perfect  interface  bond  will 
be  assumed,  this  being  the  condition  of  predominant  interest.  These  inter¬ 
face  continuity  conditions  are  automatically  satisfied  under  a  standard 
finite  element  formulation  and  need  not  be  considered  explicitly. 

In  the  usual  transverse  normal  loading  problem  (and  loading  problems 
in  general) ,  it  is  desirable  to  be  able  to  specify  the  magnitude  of  the 
applied  normal  loading  (applied  stress)  increments,  i.e.,  increments  of 
0^  and  Gy  (see  Fig.  1  or  2).  But  the  boundary  conditions  to  be  satis¬ 
fied  at  each  boundary  node  point  for  the  problem  under  study  must  be 
expressed  in  terms  of  normal  displacements  rather  than  normal  stresses, 
as  previously  explained.  This  is  significant  because  it  complicates  the 
solution  process  and  increases  the  amount  of  computation  required. 

One  method  of  achieving  specified  applied  stress  increments  at  the 
boundaries  while  satisfying  the  normal  displacement  conditions  is  to 
solve  more  than  one  boundary-value  problem  for  each  applied  stress  in¬ 
crement.  Since  the  material  response  is  approximated  to  be  linear  within 
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each  increment,  these  individual  solutions  then  can  be  superposed  to 

obtain  the  desired  result.  The  procedure  is  basically  the  same  as  that 

(2) 

used  in  previous  linear  elastic  analyses.  To  permit  the  specifica¬ 

tion  of  arbitrary  applied  stress  increments  Aa^  and  A  two  separate 
boundary-value  problems  must  be  solved  for  each  increment.  They  are 
(referring  to  Fig,  2) : 

Problem  1 

At  =  0  along  all  four  rectangular  boundaries 

Au  =  0  along  x  =  0  (boundary  node  points  remain  on  the 

coordinate  axis  because  of  symmetry) 

Au  =  1  along  x  =  a  (arbitrarily  specified  uniform  dis¬ 
placement  increment) 

Av  =  0  along  y  =  0  (boundary  node  points  remain  on  the 
coordinate  axis  because  of  symmetry) 

Av  =  0  along  y  =  b  (any  uniform  displacement  increment 
is  admissible) 

Problem  2 

Same  as  problem  1  except: 

Au  =  0  along  x  -  a 

Av  =  1  along  y  =  b 

For  each  problem,  having  solved  for  the  displacement  increment  fields 
Au1  and  Av1  (Eq.  (51)),  where  i  =  1,  2  represent  the  results  of  Problems 
1  and  2,  respectively,  the  unknown  normal  force  increment  components  at 
each  node  point  along  the  boundaries  x  =  a  and  y  =  b  can  be  computed 
from  Eq.  (50).  The  sums  of  these  force  increment  components  along  each 
boundary,  divided  by  the  corresponding  boundary  length  b  or  a,  represent 
the  average  applied  stress  increments  Acf^1  and  Aa  1  associated  with  each 
problem. 

Problems  1  and  2  then  can  be  combined  as  follows  to  obtain  the 
solution  for  arbitrarily  specified  applied  stress  increments  A and 
Act  .  Multiply  the  results  of  Problem  2  by  a  constant  a  such  that  the 
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ratio  K  of  the  desired  stress  increments  is  obtained.  (K  is  permitted 
to  vary  from  increment  to  increment  in  this  analysis.) 


-  1  -  2 

Aa  A  a  +  aAa 

k  =  — =  —i - jl. 

-  -  1  -  2 

A a  Aa  +  aAa 

XX  X 


Solving  for  a: 


KAa  1  -  Aa  1 
x  y 

-  2  -  2 

Aa  -  KAa 
y  x 


(52) 


(53) 


The  sum  of  the  Problem  1  results  plus  the  Problem  2  results  multiplied 

by  a  then  can  be  multiplied  by  a  constant  y  to  obtain  the  specific 

values  of  Aa  and  Aa  desired,  e.g. , 
x  y  5  &  * 


+  aAa  2) 

x 


(54) 


which  defines  y. 


Y  « 


Aa 


x 


-1  -  2 

Aa  +  a  a 


(55) 


x  x 

Specifically,  it  is  the  node  point  displacement  fields  of  Problems  1  and 
2  that  are  to  be  combined  in  the  above  ratio,  i.e., 


1  2 

Au  =  y(Au  +  aAu  ) 

(56) 

1  2 

Av  =  y(Av  +  aAv  ) 


All  other  quantities  of  interest  are  then  computed  using  these  node  point 
displacement  increments.  For  example,  the  strain  increment  components  in 
each  finite  element  can  be  obtained  from  Eq .  (35).  The  corresponding 
stress  increment  components  can  be  obtained  by  substituting  the  node  point 
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displacement  increment  components  into  Eq.  (36),  or,  if  the  strain  in¬ 
crement  components  have  already  been  computed,  by  substituting  these 
strain  increment  components  into  Eq.  (30). 

In  this  analysis,  the  components  A O  and  A 0  for  each  applied  stress 

x  y 

increment  can  be  specified  arbitrarily,  rather  than  only  in  a  fixed  ratio 
as  in  Ref.  1.  Because  of  this  nonproportional  loading,  local  unloading 
may  occur  even  if  the  applied  stresses  are  increased  monotonically ;  and, 
of  course,  widespread  element  unloading  will  occur  if  the  applied  stress 
increments  are  reversed  in  sign  from  previous  increments.  Thus,  it  is 
necessary  to  provide  for  unloading  response. 

The  plastic  component  of  the  total  strain  is  assumed  to  be  non- 
recoverable  and  hence  material  unloading  will  be  linearly  elastic.  Com¬ 
putationally,  the  procedure  is  as  follows.  At  the  end  of  each  load 
increment,  the  computed  octahedral  shear  stress  in  each  element  is 
compared  with  the  corresponding  value  from  the  previous  load  increment. 

If  has  decreased  in  magnitude,  the  element  is  unloading.  The  element 
material  properties  are  set  equal  to  the  elastic  values  for  the  next  and 
subsequent  load  increments,  and  held  constant  until  the  computed  value 
of  at  the  end  of  any  load  increment  again  equals  or  exceeds  the  value 
computed  for  the  increment  in  which  unloading  was  detected.  If  this 
occurs,  the  element  is  again  assigned  an  appropriate  value  of  the  elasto- 
plastic  material  parameter  2M^,  (as  defined  in  Eq .  (12))  for  the  next  load 
increment.  Since  the  load  increments  are  assumed  to  be  small  and  the 
element  stress-strain  response  is  constrained  to  follow  closely  the 
specified  piecewise  linear  curve,  errors  due  to  overshoot  and  similar 
response  lags  are  negligible. 

The  above  procedure  for  combining  Problems  1  and  2  is  used  for  all 
applied  stress  increments  except  the  first,  up  to  first  failure.  Since 
the  initial  material  response  is  linearly  elastic,  an  applied  stress  in¬ 
crement  of  arbitrary  magnitude  is  permissible  within  this  region.  Thus, 
rather  than  solve  for  a  number  of  small  stress  increments  within  the 
elastic  region,  one  increment  just  large  enough  to  stress  the  material 
of  any  one  of  the  elements  to  its  elastic  limit  value  is  computed  as  part 
of  the  solution  for  the  first  increment.  The  procedure  is  as  follows. 
Problems  1  and  2  are  solved  as  outlined  above,  and  a  value  of  a  is  computed 


(Eq.  (53))  for  the  specified  ratio  K  (Eq.  (52))  for  the  first  (elastic 
limit)  load  increment.  The  displacement  fields  of  Problems  1  and  2  are 
combined  using  Eq .  (56)  with  y  deleted  (i.e.,  as  if  y  =  1).  The  corre¬ 
sponding  element  stresses  are  then  computed  as  outlined  previously.  The 
octahedral  shear  stress  in  each  element  can  now  be  computed  from 

Eq .  (7).  The  material  of  each  element  has  an  elastic  limit  value  of 

£ 

octahedral  shear  stress  T  defined  by  its  stress-strain  curve.  Thus, 

£  0 

the  ratio  T0/T0  can  be  computed  for  each  element  and  that  element  having 
the  highest  ratio  identified.  Note  that,  in  general,  this  highest  ratio 
may  be  greater  or  less  than  unity  depending  upon  the  magnitude  of  the 
assumed  boundary  displacement  increments  used  in  Problems  1  and  2  (indi¬ 
cated  as  unity  here).  The  displacement  increment  fields  from  Eq.  (56) 
(with  y  =  1)  are  divided  by  this  ratio  to  give  the  displacement  incre¬ 
ments  corresponding  to  the  elastic  limit  loading  of  the  composite.  The 


applied  stress  increments  Ao^  and  Acr^  are  then  computed  from  these  dis¬ 
placement  increments  as  previously  outlined,  as  are  the  other  quantities 


of  interest. 


Thus,  to  summarize,  for  the  first  applied  stress  increment,  only 
the  desired  ratio  K  of  the  applied  stress  components  is  specified;  the 


specific  values  of  Ao^  and  Aa^  corresponding  to  the  elastic  limit 
computed  as  part  of  the  solution  process.  For  subsequent  applied 


are 

stress 


increments  up  to  first  failure,  values  of  Ao^  and 
solved  for. 


are  specified  and 


A  slightly  different  procedure  was  used  in  the  earlier  investigation 
described  in  Ref.  1.  The  intent  there  was  to  avoid  the  need  for  solving 
two  boundary-value  problems  for  each  loading  increment  beyond  the  first. 
The  solution  of  these  individual  boundary-value  problems  constitutes  the 
dominant  part  of  the  total  computational  time  required  (whether  it  be  the 
two  matrix  inversions  with  the  tangent  modulus  method,  or  the  two  iter¬ 
ation  processes  with  the  method  of  initial  strains).  Thus,  the  total 
computational  time  required  can  be  cut  almost  in  half  if  only  one 
boundary-value  problem  per  increment  is  necessary. 

The  first  (elastic  limit)  increment  was  obtained  as  outlined  here, 

solving  for  a  specified  ratio  of  the  applied  stress  components; 

K  =  Aa  /Aa  .  However,  unlike  this  analysis,  in  which  Aa  and  Aa  can 
y  x  .  x  y 
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be  arbitrarily  specified  positive  or  negative  values  for  each  loading 

increment,  the  ratio  K  was  held  constant  for  all  loading  increments, 

i.e.,  proportional  loading  was  assumed.  Thus,  only  a  value  of  Acr^  needed 

to  be  specified  for  each  loading  increment.  This  proportional  loading  and 

an  additional  assumption  of  monotonically  increasing  applied  stresses 

(positive  values  of  Acr^  only)  insured  that  no  local  material  unloading 

occurred.  On  these  bases,  a  predictor-corrector  technique  was  developed 

in  which  the  boundary  displacement  increment  Av  was  estimated  in  terms  of 

the  assumed  value  of  Au  and  the  results  of  the  previous  load  increment 

(thus  providing  the  reason  for  solving  for  the  first  load  increment  in  the 

usual  manner) .  The  ratio  of  the  corresponding  average  boundary  stress 

increments  Ac  /Ac  will  differ,  in  general,  from  the  desired  ratio  K,  in- 
y  x 

dicating  the  error  in  the  approximation  of  Av.  By  making  a  suitable 
correction  in  the  estimated  Av  of  the  next  increment,  however,  this  error 
was  kept  negligibly  small.  Full  details  of  the  procedure  and  the  equa¬ 
tions  used  are  presented  in  Ref.  1. 

This  procedure  worked  well,  being  adequately  stable  in  loading  the 
composite  to  first  failure.  However,  several  subsequent  attempts  to 
adopt  a  similar  procedure  for  post-first-failure  applied  stress  increments 
and  the  related  adjustment  increments  (which  will  be  described  in  the 
next  subsection)  were  unsuccessful.  Attempts  to  correct  for  even  small 
errors  in  Av  resulted  in  serious  instabilities.  The  amount  of  computa¬ 
tion  time  required  and  the  possible  errors  introduced  in  attempting  to 
switch  from  a  predictor-corrector  method  before  first  failure  to  a  two- 
problem  solution  after  first  failure  led  to  the  abandonment  of  this 
approach.  The  slight  savings  in  computation  time  were  more  than  offset 
by  the  errors  and  uncertainties  introduced.  Thus,  a  two-problem  incre¬ 
mental  solution  has  been  adopted  in  the  current  investigation,  for  all 
applied  stress  (and  adjustment)  increments. 

At  some  point  as  the  incremental  loading  procedure  is  continued,  the 
octahedral  shear  stress  and/or  the  octahedral  plastic  shear  strain  in  one 
or  more  elements  will  exceed  the  ultimate  values  and  f°r  that 

element  as  defined  by  the  material  stress-strain  curve.  Note  that  in 
theory  the  material  response  should  follow  the  stress-strain  curve 
exactly,  and  thus  TU  and  e  would  be  reached  simultaneously.  When 
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approximating  the  curve  by  linear  segments  as  is  done  here,  however,  it 

is  possible  that  one  limiting  value  may  be  exceeded  without  exceeding 

the  other.  The  criterion  for  element  failure  used  here  will  be  that 

both  TU  and  £  must  be  reached, 

o  o  ,  v 

Tf  *rU  ^  are  exceeded  in  at  least  one  element  for  a 


If  both  T  and  £ 
o  o 


given  applied  stress  increment,  the  specified  values  of  A and  Ac^  for 

that  increment  are  then  reduced  to  values  just  sufficient  to  cause  fail- 

ure  of  only  one  element.  If  T°  or  £  (but  not  both)  is  exceeded  in 

oo 

at  least  one  element,  the  specified  values  of  Aa  and  Ac  for  that  in- 

»  r  x  y 

.  u 

crement  are  then  reduced  to  values  just  sufficient  to  cause  either  T 

o 

or  £  to  be  reached  in  just  one  element.  In  this  latter  case,  an 

o 

element  failure  has  not  occurred,  however,  and  the  next  load  increment 

is  applied  and  the  process  repeated. 

This  reduction  of  the  applied  stress  increments  Ao^  and  Ac^  for 

load  increments  beyond  the  first  (linearly  elastic)  increment  cannot  be 

achieved  by  a  simple  linear  adjustment  such  as  was  used  in  obtaining  the 

elastic  limit  values  in  the  first  load  increment.  For  an  applied  stress 

increment,  both  T  and  £  ^  (Eqs .  (7)  and  (6),  respectively)  are  func- 
°  0  2 

tions  of  terms  of  the  form  (a  +  Ac  )  ,  where  a  is  a  stress  component 

XX  x 

in  an  element  at  the  beginning  of  the  increment  and  Ac  represents  the 

o  o 

incremental  increase.  Thus,  terms  of  the  form  ax  +  2cxAcx  +  Acx  must 
be  scaled  down  to  achieve  the  desired  value  of  or  £0^P^U.  That  is, 
a  quadratic  equation  must  be  solved  in  determining  the  scale  factor  in 
each  case. 

These  quadratic  equations  reduce  to  a  single  linear  relation  for 
the  first  load  increment  since  the  stress  components  at  the  beginning 
of  the  increment  (those  represented  by  in  the  above  expression)  are 
zero,  making  the  previously  discussed  linear  scaling  procedure  (Eqs. 

(54)  and  (55))  valid. 

Equaling  or  exceeding  both  and  f°r  one  element  defines 

first  failure  of  the  composite  body.  This  local  failure  must  then  be 
accounted  for,  as  will  be  discussed  in  the  next  subsection. 
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CRACK  INITIATION  AND  PROPAGATION  TO  TOTAL  FAILURE 

When  the  material  in  a  local  region  is  stressed  to  its  ultimate 
value,  as  defined  by  the  stress-strain  response  for  that  constituent 
material,  it  will  fail.  If  this  occurs  in  a  region  of  high  stress 
gradients  such  as  typically  exist  in  a  composite  material,  the  surround¬ 
ing  material  may  be  able  to  absorb  the  redistribution  of  stresses  caused 
by  the  local  failure  without  additional  failure  occurring  at  that  level 
of  applied  stress.  If  this  is  the  case,  the  initial  failure  becomes  a 
local  discontinuity--a  crack — within  the  material.  As  additional  load¬ 
ing  is  applied,  this  crack  may  grow  in  size  (propagate),  or  additional 
cracks  may  be  initiated  in  other  regions  of  high  stress  within  the 
material.  Eventually,  the  remaining  material  will  be  unable  to  carry 
the  applied  loads  and  total  failure — a  crack  completely  across  the 
material  cross  section — will  result. 

Thus,  it  is  first  necessary  to  model  the  crack.  Since  the  material 
region  is  represented  by  elements  of  finite  size,  in  each  of  which  the 
stresses  and  strains  are  assumed  constant,  the  imminence  of  a  local 
failure  is  identified  with  a  particular  finite  element.  An  obvious 
method,  and  that  used  here,  is  to  subsequently  set  the  material  proper¬ 
ties  of  that  element  equal  to  zero,  i.e.,  the  stresses  and  the  stiffness 
properties  within  the  element  are  reduced  to  zero.  Thus,  the  crack  is 
assumed  to  have  the  dimensions  of  the  failed  element.  In  an  actual 
material,  the  width  of  the  crack  will  typically  be  considerably  smaller. 

This  "failed  element"  approximation  of  a  crack  has  two  implications. 
A  finite  amount  of  material  is  assumed  to  be  removed  from  the  system, 
which  is  not  actually  the  case ,  and  the  crack  is  not  likely  to  completely 
close  up  under  subsequent  loadings  because  of  its  exaggerated  width. 
Should  an  actual  crack  subsequently  close  up,  it  could  support  compres¬ 
sive  normal  forces  and  shear  forces  across  its  surfaces.  Of  course,  as 
the  finite  elements  are  made  smaller  and  smaller,  the  model  representa¬ 
tion  of  an  actual  crack  improves. 

There  are  alternative  methods  of  modeling  the  crack  that  do  not 
involve  the  undesirable  features  noted  above;  one  of  them  is  outlined 
in  Appendix  B.  However,  practical  difficulties  in  implementing  such 
methods  have  thus  far  prevented  their  use.  In  this  investigation,  we 
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use  the  "failed  element"  concept  of  setting  the  element  material  proper¬ 
ties  equal  to  zero  and  defining  the  resulting  void  as  the  crack.  As 
will  be  demonstrated  in  Secs.  IV  and  V,  this  method  has  proven  to  be 
quite  satisfactory. 

The  finite  element  method  involves  the  maintenance  of  force  equi¬ 
librium  at  each  node  point  in  the  array.  Each  element  surrounding  a 
given  node  point  exerts  a  force  on  that  node,  and  the  magnitude  of  these 
forces  is  proportional  to  the  stresses  in  the  elements,  as  given  by 
Eq .  (47).  Thus,  to  represent  the  unloading  due  to  failure  of  an  ele¬ 
ment,  equal  and  opposite  force  components  are  applied  at  the  correspond¬ 
ing  three  node  points.  This  reduces  the  node  point  force  contribution 
of  that  element  to  zero;  since  the  material  stiffness  properties  of  the 
element  are  simultaneously  set  equal  to  zero,  no  stresses  will  be  de¬ 
veloped  within  the  element  during  subsequent  applied  stress  increments. 

These  negating  node  point  forces  will  cause  a  complete  redistribu¬ 
tion  of  displacements,  strains,  and  stresses  through  the  composite 
material.  Thus,  after  each  applied  stress  increment  in  which  an  element 
failure  is  indicated,  an  "adjustment"  increment  must  be  performed  to 
apply  these  negating  force  components.  The  effect  of  an  element  failure 
during  an  applied  stress  increment  is  hence  accounted  for  in  the  sub¬ 
sequent  adjustment  increment. 

As  in  the  case  of  the  applied  stress  increments,  certain  boundary 
conditions  must  be  satisfied  during  the  adjustment  increment.  These 
conditions  are  zero  shear  stresses  and  uniform  displacements  along  the 
rectangular  boundaries  of  the  first  quadrant  (see  Fig.  2). 

A  choice  can  be  made  at  this  point  between  maintaining  constant 
average  boundary  stresses  (A a  =  Act  =0)  during  the  adjustment  incre- 

x  y  _  _ 

ment ,  constant  boundary  displacements  (Au  =  Av  =  0) ,  or  a  combination 
of  these  conditions.  In  this  study,  a  combination  has  been  used  that 
simulates  the  conditions  of  a  typical  uniaxial  test  using  a  constant 
deformation  testing  machine.  The  test  specimen  is  held  fixed  in  the 
loading  direction  (along  the  x-axis)  during  the  adjustment  process, 
with  its  lateral  dimension  (along  the  y-direction)  free  to  expand  or 
contract  while  maintaining  a  constant  average  stress  (which,  by  defi¬ 
nition,  is  zero  for  the  special  case  of  uniaxial  loading).  This 
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combination  requires  the  solution  of  two  separate  boundary- value  prob¬ 
lems,  as  does  the  applied  stress  increments.  The  required  boundary 
conditions  are  different,  but  the  method  of  combining  them  is  basically 

the  same.  The  previously  defined  Problem  2  is  unchanged.  However, 

Problem  1  is  replaced  by  Problem  3,  defined  as  follows: 

Problem  3 

At  =0  along  all  four  rectangular  boundaries 
Au  =  0  along  x  =  0  and  x  =  a 

Av  =  0  along  y  =  0  and  y  =  b 

Force  components  of  equal  magnitude  but  opposite  sign  from  those 

given  by  Eq .  (47)  are  applied  at  the  node  points  of  the  failed 

element . 

Problems  2  and  3  are  then  combined  to  obtain  the  following  adjust¬ 
ment  increment  results: 

Au  =  0  along  x  =  a  (obtained  automatically  since  Au  =  0 
for  both  Problems  2  and  3) 

AOy  =  0  along  y  =  b 

The  displacement,  strain,  and  stress  increment  fields  are  then  added 
to  the  previously  accumulated  sums  in  the  usual  manner. 

The  redistributions  of  stresses  due  to  this  adjustment  increment 
may  cause  one  or  more  additional  elements  to  exceed  their  ultimate 
octahedral  shear  stress  and  octahedral  plastic  shear  strain  values, 
i.e.,  to  fail.  If  this  happens,  an  additional  adjustment  increment 
is  solved  for.  This  process  is  repeated  until  no  additional  elements 
fail  during  an  adjustment  increment.  The  next  applied  stress  increment 
is  then  applied  in  the  usual  manner. 

As  previously  noted,  adjustment  increment  boundary  conditions  other 
than  those  selected  for  this  analysis  can  also  be  used,  the  choice  de¬ 
pending  upon  the  physical  problem  to  be  simulated. 

For  example,  to  simulate  a  constant  deformation  (e.g.,  screw-driven) 
biaxial  loading  testing  machine,  the  adjustment  increment  boundary 
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conditions  Au  =  0  along  x  -  a  and  Av  =  0  along  y  =  b  would  be  appro¬ 
priate.  Since  this  corresponds  to  Problem  3  outlined  above,  only  one 
boundary-value  problem  would  have  to  be  solved  per  adjustment  increment, 
and  no  combining  of  individual  problems  is  necessary. 

To  simulate  a  constant  load  (e.g.,  hydraulic  or  dead  load)  biaxial 
loading  testing  machine,  the  adjustment  increment  boundary  conditions 

Aa  =0  along  x  =  a  and  Ao  =0  along  y  =  b  would  be  appropriate.  This 
x  y 

would  require  the  combination  of  solutions  for  Problems  1,  2,  and  3  for 
each  adjustment  increment.  A  procedure  for  combining  three  problems 
is  discussed  in  detail  in  Ref.  12  (in  considering  the  influence  of  a 
uniform  temperature  change)  and  is  not  repeated  here. 

These  and  similar  other  adjustment  increment  boundary  conditions 
can  be  alternatively  readily  incorporated  as  desired. 

Total  failure  of  the  composite  occurs  when  a  continuous  band  of 
failed  elements  extends  completely  across  the  first  quadrant,  either 
from  the  x  =  0  boundary  to  the  x  =  a  boundary,  or  from  the  y  =  0  boun¬ 
dary  to  the  y  =  b  boundary.  This  necessarily  corresponds  to  either 

O  =  0  or  a  =  0,  with  no  further  change  in  O  or  O  possible.  In  sum- 
y  x  &  y  x  r 

mary ,  a  crack  has  propagated  completely  across  the  composite  material, 
causing  total  failure. 
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III.  COMPUTATIONAL  ASPECTS 


This  section  is  intended  for  those  who  are  interested  in  the  gen¬ 
eral  characteristics  of  the  numerical  solution  procedure. 

Constant  strain  triangular  elements,  as  defined  in  Sec.  II,  are  used 
to  depict  the  region  of  interest.  The  computer  program  is  presently 
dimensioned  to  allow  the  use  of  up  to  350  elements  and  200  node  points. 
The  size  of  the  element  array,  and  specifically  the  number  of  node  points, 
directly  influences  the  computer  time  required  per  solution  increment 
since  the  solution  is  obtained  in  terms  of  node  point  displacement  com¬ 
ponents.  Thus,  a  trade-off  must  always  be  made  between  improving  the 
resolution  by  using  a  larger  number  of  elements  and  node  points,  and 
keeping  the  required  computation  time  within  acceptable  bounds.  The 
limits  on  the  array  size  have  been  more  than  adequate  to  date.  Much  can 
be  done  to  improve  the  resolution,  without  increasing  computation  time, 
by  carefully  constructing  the  finite  element  grid. 

The  composite  material  configuration  that  can  be  portrayed  within 
the  region  being  analyzed  (the  first  quadrant  of  the  typical  repeating 
unit  indicated  in  Fig.  2)  is  completely  arbitrary.  Any  filament  shape 
that  can  be  represented  by  a  group  of  triangular  elements  is  admissible, 
and  more  than  one  filament  can  be  included  in  the  first  quadrant.  A 
void  can  be  expressed  by  setting  the  Young's  modulus  equal  to  zero  for 
those  elements  that  represent  its  shape.  The  program  is  dimensioned 
to  handle  as  many  as  five  different  materials. 

Since  an  elastoplastic  analysis  is  being  performed,  the  complete 
stress— s train  response  to  total  failure  for  each  constituent  material 
must  be  given.  The  analysis  permits  any  monotonic  response  (i.e.  , 
single-valued  relation  between  stress  and  strain)  to  be  specified.  The 
stress-strain  curve  for  each  material  is  input  point  by  point,  usually 
using  actual  experimental  data.  The  applied  stress  components  Aa^  and 
Act  can  be  arbitrarily  specified  for  each  solution  increment.  The  pro¬ 
gram  is  dimensioned  to  handle  as  many  as  50  pre-first-failure  load 
increments  and  50  post-first-failure  load  increments.  A  typical  prob¬ 
lem  may  require  10  to  12  pre-first-failure  load  increments  and  12  to 
15  post-first-failure  load  increments. 
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The  finite  element  computer  program  provides  a  printout  of  results 
for  each  applied  stress  increment  and  each  adjustment  increment,  namely, 
the  increments  of  average  composite  stress  actually  solved  for,  the 
corresponding  increments  of  composite  strain,  the  accumulated  sums  of 
these  quantities,  and  the  accumulated  node  point  displacement  compon¬ 


ents.  It  also  prints  out  for  each  element  the  accumulated  values  of 

the  strains  £  ,  £  ,  y  ,  z  ^  and  the  following  normalized  quantities: 
x  y  xy  o 

the  stress  components  O  /d  y  a  /a  9  o  /a  9  x  /a  ;  the  octahedral  shear 
^  xx’  yx  zx  xyx 

stress  T  /t  ;  the  maximum  principal  stress  O^/d  ;  and  the  maximum  prin- 
oo  lx 

cipal  strain  £-,/£  .  Furthermore,  it  prints  out  the  numbers  of  the  ele- 

ments  for  which  either  TU  or  £  P  U  has  been  exceeded  during  the  solution 

o  o 

increment . 


In  addition  to  this  printed  output,  the  program  prepares  a  tape 
containing  those  results  that  are  required  for  input  to  the  contour 
plotting  computer  program,  used  by  a  Stromberg-DatagraphiX  S-D  4060 
plotter . 


The  finite  element  grid,  which  is  input  node  by  node  into  the 
finite  element  computer  program,  can  be  plotted  with  and  without  ele¬ 
ment  numbers  being  printed.  This  serves  as  a  ready  visual  check  for 
the  correctness  of  the  input  data.  Contour  plots  of  any  of  the  follow¬ 
ing  element  quantities  can  be  specified  (normalized  as  indicated) : 
maximum  principal  stress  O./d  ,  maximum  principal  strain  £- /e  ,  octa- 

^  \  l  X  g 

hedral  plastic  shear  strain  £  ^  ,  and  octahedral  shear  stress  T  /x  . 

o  o  o 

Values  of  the  contours  to  be  plotted  are  specified  as  input  data. 

Failed  elements  are  also  indicated  on  the  contour  plots.  These 
elements  are  outlined  and  shaded,  making  it  easy  to  follow  visually  the 
progress  of  a  crack  on  plots  of  successive  solution  increments. 
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IV.  NUMERICAL  RESULTS 


Because  the  primary  goal  of  this  investigation  has  been  to  develop 
a  method  of  analysis  and  to  write  an  associated  computer  program,  most 
of  the  concentration  has  been  on  attempting  to  improve  computational 
techniques  and  solution  accuracy.  Only  very  limited  numerical  results 
have  been  „ obtained  to  date. 

However,  in  the  process  of  refining  the  detailed  aspects  of  the 
numerical  solution  procedure  and  correcting  programming  errors,  a  sig¬ 
nificant  number  of  computer  runs  have  been  made.  This  development 
effort  involved  the  study  of  effects  of  changing  finite  element  grid 
patterns  in  local  regions,  modifying  boundary  conditions,  reshaping 
the  governing  stress-strain  curves  of  the  constituent  materials,  vary¬ 
ing  the  size  of  the  applied  stress  increments,  and  so  forth.  Thus,  the 
few  examples  presented  here  and  the  numerical  analysis  used  to  obtain 
their  solutions  constitute  a  summary  of  the  knowledge  gained  to  date. 

There  was  no  single  computer  run  that  was  completely  satisfactory 
in  terms  of  modeling  all  aspects  of  the  composite  material  response 
from  the  first  (elastic  limit)  applied  stress  increment,  to  first  fail¬ 
ure,  and  then  to  complete  crack  propagation  and  total  composite  failure. 
Each  of  these  separate  phases,  however,  was  successfully  run  at  various 
times  during  the  computer  program  development  process.  An  inherent 
numerical  inaccuracy  of  the  constant  strain  triangular  finite  element 
technique,  which  remains  to  be  resolved  (as  will  be  discussed  in  the 
next  section) ,  combined  with  a  lack  of  available  computer  time  has 
limited  the  number  of  numerical  results  obtained. 

The  examples  presented  here  are  intended  to  demonstrate  the  type 
of  information  which  can  be  obtained,  even  though  the  results  are  too 
limited  to  provide  conclusive  predictions  of  actual  crack  propagation 
and  composite  failure.  Many  more  computer  runs  will  be  required  to 
achieve  this. 

A  boron-aluminum  composite  was  assumed  in  most  of  the  numerical 
examples  considered.  The  distinct  change  in  slope  of  the  stress-strain 
curve  beyond  the  elastic  limit  for  the  annealed  6061  aluminum  alloy 
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matrix  material  that  was  selected  offered  a  good  test  of  the  elasto- 

plastic  response  representation  of  the  analysis.  Actual  experimental 

* 

data  provided  by  TRW-Cleveland  were  used;  the  stress-strain  curve  ob¬ 
tained  from  a  uniaxial  tensile  test  of  this  matrix  material  is  shown 

in  Fig.  6.  The  complete  curve  is  input  point  by  point  into  the  computer 

6  2 

program.  The  measured  Young's  modulus  was  9.4  x  10  lb/in.  ,  the  elastic 

2 

limit  stress  was  6280  lb/in.  ,  and  failure  occurred  at  an  ultimate  stress 
2 

of  18,900  lb/in.  and  a  strain  of  30  percent.  The  boron  filaments  were 

assumed  to  exhibit  linearly  elastic  response  to  failure,  based  upon 

longitudinal  tensile  tests  of  individual  filaments.  Their  measured 

6  2 

Young's  modulus  was  55.0  x  10  Ib/in.  ,  and  an  ultimate  strength  of 
2 

500,000  lb/in.  was  assumed  as  being  sufficiently  high  to  insure  that 
failure  would  occur  in  the  matrix  only.  This  was  the  observed  failure 
mode  in  the  actual  composite  material  being  modeled.  The  Poisson's 
ratios  for  the  aluminum  and  boron  were  assumed  to  be  0.32  and  0.20, 
respectively . 

Complete  stress-strain  curves  to  failure  for  a  number  of  uniaxial, 

transverse  normal  tensile  tests  of  composites  using  these  constituent 

materials  were  provided  by  TRW-Cleveland,  in  the  form  of  strip  chart 
& 

records.  Photomicrographs  of  the  test  specimens  before  and  after 

(13) 

failure  were  also  available.  All  specimens  using  the  annealed  alum¬ 

inum  matrix  failed  in  the  matrix.  Other  composite  specimens,  using  the 
same  6061  aluminum  alloy  but  in  the  heat-treated  condition,  also  exhibited 
failures  in  the  filaments.  The  alloy  is  much  stronger  in  the  heat-treated 
condition,  as  indicated  by  the  experimental  curve  in  Fig.  6.  A  complete 
discussion  of  the  test  procedures,  specimen  designs,  material  processing, 
and  experimental  program  in  general  may  be  found  in  Ref.  13. 

One  specific  test  result,  TRW  Specimen  4Fq2/2,  that  was  selected  to 
be  modeled  by  the  analysis  is  an  annealed  specimen  containing  40  volume 
percent  (v^  =  40  percent)  of  boron  filaments.  This  specimen  was  selected 
because  it  exhibited  a  stress-strain  curve  that  represented  well  the 
data  available.  The  filament  packing  was  very  regular  and  in  approx¬ 
imately  a  square  array.  Photomicrographs  are  included  in  Ref.  13. 


Private  communication  of  unpublished  experimental  data  from  I.  J. 
Toth,  TRW,  Inc.,  Cleveland,  Ohio. 
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Fig.  6 —  Uniaxial  tensile  stress-strain  curves,  6061  aluminum  alloy 
matrix  materials  (data  from  Ref.  13) 
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The  finite  element  grid  shown  in  Fig.  3  was  constructed  to  model 
this  geometry,  boron  filaments  being  circular  in  cross  section.  This 
is  a  relatively  large  grid,  containing  176  nodes  and  304  elements,  and 
it  models  the  matrix  region  very  well.  A  coarse  element  grid  was  used 
in  the  filament  region  since  the  boron  filaments  remain  elastic  and  do 
not  fail.  This  conserved  node  points  and  should  not  significantly  af¬ 
fect  solution  accuracy. 

The  actual  composite  stress-strain  response  obtained  experimentally 

by  TRW  is  indicated  by  the  solid  curve  in  Fig.  7.  This  was  a  uniaxial, 

transverse  normal  tensile  test,  corresponding  to  a  a  component  of 

x 

applied  stress,  as  indicated  in  Fig.  1 ;  is  zero.  One  set  of  numer¬ 
ical  results  of  the  present  analysis,  computer  run  8A-10-17,  is  illus¬ 
trated  by  the  dashed  curve  in  Fig.  7. 

First  yielding  was  predicted  to  occur  in  the  matrix  along  the 
x-axis  midway  between  adjacent  filaments — in  element  79  in  the  lower 
right  corner  of  the  grid  of  Fig,  3.  This  is  consistent  with  the  re¬ 
sults  obtained  in  Ref.  1.  As  noted  in  Fig.  7,  first  yielding  was  pre- 

2 

dieted  to  occur  at  an  applied  stress  of  5539  lb/in.  and  a  composite 

strain  of  0.0296  percent,  resulting  in  a  predicted  composite  elastic 

6  2 

modulus  of  18.7  x  10  lb/in.  .  The  experimentally  measured  elastic 

6  2 

modulus  was  17.3  x  10  lb/in.  .  As  indicated  in  Fig.  7,  the  initial 

portions  of  the  two  curves  are  in  good  agreement. 

2 

The  composite  was  then  loaded  in  500  lb/in.  increments  up  to  first 

failure,  which  was  predicted  to  occur  in  element  139,  an  element  of 

matrix  material  located  at  the  filament-matrix  interface  (see  Fig.  3). 

This  occurred  during  applied  stress  increment  no.  29,  at  an  applied 

2 

stress  of  19,519  lb/in.  .  This  point  is  marked  first  failure  in  Fig.  7. 

The  failed  element  is  outlined  in  Fig.  8;  the  number  inside  the  element 

indicates  the  applied  stress  increment  during  which  failure  occurred. 

Octahedral  shear  stress  contours  (which  have  been  normalized  by  dividing 

by  the  value  of  the  octahedral  shear  stress  at  which  yielding  of  the 

£  2 

matrix  material  occurs,  i.e.,  T  =  2960  lb/in.  )  are  also  indicated. 

0  2 
Matrix  failure  occurs  at  an  octahedral  shear  stress  of  8910  Ib/in.  , 

u  £ 

corresponding  to  a  normalized  value  of  t  /t  =  8910/2960  =  3.01.  Thus. 

o  o  ’ 

the  contours  having  the  value  1.0  indicate  the  initiation  of  yielding. 
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Fig  .  7—  Predicted  and  experimental  stress-strain  curves 
Computer  run  8A-  10-17 
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Fig  .  8 —  First  element  failure  and  octahedral  shear  stress  contours  (normalized) 
Computer  run  8A-10-17,  applied  stress  increment  29 
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and  those  having  the  value  2.0  indicate  a  stress  level  approximately 
midway  between  yielding  and  failure. 

The  first  element  failure,  as  indicated  in  Fig.  8,  resulted  in  a 
decrease  of  196  lb/in.  in  the  applied  stress  during  the  immediately 
subsequent  adjustment  increment  (which  is  always  required  when  an  ele¬ 
ment  fails) ,  and  this  adjustment  did  not  cause  additional  elements  to 
fail . 

2 

Beyond  first  failure,  200  lb/in.  applied  stress  increments  were 
used.  During  applied  stress  increment  no.  30,  the  plastic  octahedral 

shear  strain  ultimate  £  (but  not  the  octahedral  shear  stress  ulti- 

°  2 
mate  TU)  was  exceeded  in  element  81  (see  Fig.  3).  This  200  lb/in. 
o 

applied  stress  increment  was  thus  automatically  normalized  to  give  the 

value  £  for  element  81,  which  corresponded  to  an  actual  applied 

0  .  2 
stress  increment  of  87  lb/in.  . 

During  the  next  applied  stress  increment,  element  81  failed,  at 

2 

a  stress  increment  of  only  24  lb/in.  .  The  subsequent  adjustment  in¬ 
crement  caused  a  decrease  in  applied  stress  of  85  lb/in.  ,  and  (but 

not  £  was  exceeded  in  elements  78  and  79.  These  two  elements 

o 

failed  during  the  next  applied  stress  increment,  increment  no.  32, 

When  this  occurred,  a  total  of  four  adjustment  increments  were  required 
to  again  achieve  equilibrium.  During  the  first  adjustment  increment, 
element  93  failed.  Thus,  a  second  adjustment  increment  was  required, 
during  which  element  82  failed.  This  required  a  third  adjustment,  dur¬ 
ing  which  elements  92,  94,  and  105  failed.  A  stable  condition  was 
achieved  during  the  fourth  adjustment  increment,  although  had  been 
exceeded  in  element  117.  At  this  point  the  applied  stress  had  decreased 
to  17,282  lb / In . 2 . 

During  the  next  applied  stress  increment,  increment  no.  33,  element 

2  2 

117  failed  after  only  10  lb/in.  of  the  allowable  200  lb/in.  increment 

had  been  applied.  This  in  turn  caused  element  106  to  fail  during  the 

subsequent  adjustment  increment,  requiring  a  second  adjustment  increment. 

2 

This  restored  equilibrium,  at  an  applied  stress  of  16,246  lb/in.  »  and 
completed  the  first  major  decrease  in  applied  stress,  as  indicated  in 
Fig.  7.  The  elements  which  had  failed  at  this  point  are  indicated  in 
Fig.  9.  When  the  number  within  an  element  is  hyphenated,  the  first 
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Fig. 9 —  Failed  elements  and  octahedral  shear  stress  contours  ( normal ized ) 
Computer  run  8A-10-17,  applied  stress  increment  33-2 
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digits  indicate  the  applied  stress  increment  number,  and  those  after 
the  hyphen  indicate  the  adjustment  increment  within  that  applied  stress 
increment  during  which  the  element  failed.  When  the  number  is  not 
hyphenated,  it  indicates  that  that  element  failed  during  the  applied 
stress  increment  itself.  Shading  identifies  the  element  that  failed 
in  the  previous  figure. 

Seven  additional  200  lb/in.  applied  stress  increments,  nos.  34 

through  4CL,  were  then  employed;  no  additional  elements  failed.  At 
2 

83  lb/in.  of  increment  no.  41,  element  127  failed.  This  in  turn  caused 

element  118  to  fail,  and  equilibrium  was  restored  at  an  applied  stress 
2 

of  16,707  lb/in.  .  This  decrease  is  represented  by  the  second  sharp 
drop  in  applied  stress  in  Fig.  7,  at  a  composite  strain  of  0.258  percent. 

The  process  of  applying  stress  increments  and  causing  additional 
elements  to  fail  was  continued  as  outlined  above,  resulting  in  the  jagged 
stress-strain  response  indicated  in  Fig.  7.  During  this  time,  the  crack 
initially  continued  to  propagate  upward  along  the  right  boundary  of  the 
first  quadrant,  and  then  a  second  crack  began  to  propagate  along  the 
interface.  At  applied  stress  increment  no.  72,  when  element  208  failed 
(see  Fig.  3),  the  applied  stress  dropped  suddenly.  This  point  is  noted 
in  Fig.  7.  Elements  which  had  failed  at  this  time  are  indicated  in 
Fig.  10. 

Element  207  failed  immediately  during  the  next  applied  stress  incre¬ 
ment,  marking  the  point  at  which  the  crack  turned  away  from  the  interface, 
propagating  out  into  the  matrix  material  region.  During  the  first  nine 
adjustment  increments  following  this  applied  stress  increment,  the  crack 
continued  to  propagate  upward  and  at  the  same  time  it  connected  with  the 
original  crack,  as  indicated  in  Fig.  11.  At  this  point  a  numerical  in¬ 
stability  interrupted  the  solution  process. 

In  view  of  later  results,  it  is  fairly  obvious  that  the  crack  ex¬ 
tending  out  into  the  matrix  from  the  interface  would  have  continued  to 
propagate  upward  to  the  top  boundary  of  the  material  region  if  the  solu¬ 
tion  process  had  been  continued.  This  would  have  constituted  total 
failure.  Looking  again  at  Fig.  7 ,  it  is  anticipated  that  the  predicted 
stress-strain  curve  would  have  continued  downward  with  only  minor  fluc¬ 
tuations.  For  all  practical  purposes,  failure  had  effectively  occurred 
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Fig.  10 —  Failed  elements  and  octahedral  shear  stress  contours  ( normal  ized ) 
Computer  run  8A-10-17,  applied  stress  increment  72 
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Fig.  11 —  Failed  elements  and  octahedral  shear  stress  contours  (normalized) 
Computer  run  8A-10-17,  applied  stress  increment  73 -V 
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during  applied  stress  increment  no.  66,  when  the  crack  began  to  propa¬ 
gate  along  the  interface.  Beyond  this  point  the  applied  stress  de¬ 
creased  rapidly,  with  only  minor  pauses. 

Although  the  first  crack  began  at  the  interface  (during  applied 
stress  increment  no.  29),  the  major  crack  propagation  was  initially 
along  the  right  boundary  (up  the  midline  between  adjacent  filaments). 
Actually,  the  fact  that  first  failure  occurred  in  element  139  may  have 
been  an  anomaly ,  caused  by  the  distorted  shape  of  element  15  in  the 
filament  directly  across  the  interface  (see  Fig.  3).  This  distortion 
was  introduced  while  minimizing  the  matrix  bandwidth  during  construc¬ 
tion  of  the  grid  array.  Element  81  was  also  very  near  to  failure  when 
element  139  failed.  At  applied  stress  increment  no.  66,  during  which 
the  second  element  near  the  interface  failed,  the  magnitude  of  the 
applied  stress  was  17,757  Ib/in.  ,  and  the  composite  strain  was  0.388 
percent.  The  composite  was  already  near  total  failure. 

The  composite  stress-strain  curve  presented  in  Fig.  7  is  exactly 
as  obtained  in  computer  run  8A-10-17 .  Later,  however,  several  relatively 
minor  errors  in  the  computer  program  were  discovered  and  corrected.  The 
principal  effects  of  these  corrections  were  to  approximately  double  the 
plastic  component  of  the  composite  strain,  and  to  decrease  the  applied 
stress  at  which  first  failure  occurred.  Because  of  a  numerical  insta¬ 
bility  problem  that  is  yet  to  be  resolved  (as  will  be  discussed  in  the 
next  section)  and  to  conserve  available  computer  time,  a  complete  rerun 
of  this  example  problem  was  not  attempted.  However,  a  computer  rerun 
(8C-1)  up  to  first  failure  was  made.  The  result  is  plotted  in  Fig.  12 
along  with  the  original  post-first-failure  results,  which  have  been 
approximately  corrected  for  the  previous  error  by  doubling  the  plastic 
component  of  the  composite  strain  values.  Thus  it  must  be  emphasized 
that  Fig.  12  does  not  represent  an  actual  result  of  a  single  computer  run, 
but  rather  is  an  estimate  of  the  type  of  behavior  that  can  reasonably  be 
expected  in  subsequent  investigations. 

The  postulated  stress-strain  curve  of  Fig.  12  is  clearly  in  good 
agreement  with  the  experimental  data.  One  fact  that  should  be  kept  in 
mind,  however,  is  that  the  experimental  data  represent  the  results  of  a 
test  of  a  composite  containing  hundreds  of  individual  filaments  that  were 
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Composite  strain  (percent) 


Pig  .  12—  Predicted  and  experimental  stress-strain  curves 
Combined  computer  runs  8  C  ■  1  and  8A-10-17 
(adjusted  strains) 
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not  all  perfectly  spaced  in  a  square  array.  It  would  be  very  informa¬ 
tive  to  determine  the  influence  of  minor  variations  of  filament  spacing 
on  both  the  predicted  composite  stress-strain  response  and  the  actual 
crack  propagation  pattern.  The  experimental  curve  represents  a  sta¬ 
tistical  average  of  such  a  variation. 

The  first  example  discussed,  computer  run  8A-10-17,  was  unique  in 
that  a  very  large  number  of  applied  stress  increments  were  required,  a 

total  of  73.  This  was  necessitated  by  the  relatively  small  size  of  the 

2 

applied  stress  increments  specified — 500  lb/in.  prior  to  first  failure, 
2 

and  200  lb/in.  after  first  failure.  A  study  of  the  nearly  50  computer 
runs  that  were  made  during  the  investigation  have  somewhat  quantified 
what  is  intuitively  obvious:  the  numerical  results  are  influenced  by 
the  size  of  the  increments  used. 

For  example,  computer  run  8C-1  used  larger  applied  stress  incre¬ 
ments  that  did  computer  run  8A-10-17.  The  pre-first- failure  applied 

2 

stress  increments  were  specified  as  500  lb/in.  for  increment  nos.  2 

2  2 
and  3,  1000  lb/in.  for  increment  no.  4,  and  1500  lb /in.  for  all  addi¬ 
tional  pre-first-failure  increments  required.  (The  magnitude  of  the 
first  (elastic  limit)  applied  stress  increment  is  solved  for  in  the 

analysis  directly,  as  previously  discussed.)  The  computer  program  was 

2 

then  immediately  rerun  (802)  using  all  500  Ib/in.  applied  stress  in¬ 
crements.  Although  it  terminated  prior  to  first  failure  because  of 
a  numerical  instability,  the  predicted  stress-strain  curves  as  obtained 
from  runs  801  and  802,  up  to  the  point  where  the  instability  occurred 
in  run  802,  were  somewhat  different.  The  applied  stress  at  a  given 
composite  strain  was  slightly  lower  when  the  larger  increments  were 
used.  Also,  it  was  indicated  that  the  applied  stress  at  which  first 
failure  would  have  occurred  in  computer  run  8C-2  would  have  been  some¬ 
what  higher  than  was  actually  obtained  in  computer  run  801  because 
local  stress  concentrations  within  the  composite  tended  to  be  lower. 
(Probable  reasons  for  this  will  be  discussed  in  conjunction  with  the 
next  example.)  It  would  probably  have  been  more  accurate  to  use  run 
8C-2  as  the  first  portion  of  the  predicted  curve  in  Fig.  12  (to  first 
failure) ,  if  a  numerical  instability  had  not  occurred  just  prior  to 
that  point,  making  the  stress  at  first  failure  indeterminate. 
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The  influence  of  applied  stress  increment  size  was  also  demonstrated 
in  another  computer  run,  8A-19-25.  It  was  run  immediately  after  8A-10-17, 
with  only  very  minor  changes  in  the  computer  program.  The  significant 
difference  between  these  two  runs  was  in  the  size  of  the  applied  stress 
increments  specified.  Computer  run  8A-19-25  used  larger  increments.  The 
pre-first-failure  applied  stress  increments  were  500  lb/in.  for  increment 

ry  2 

nos.  2  and  3,  1000  lb/in.  for  increment  no.  4,  and  1500  lb/in.  for  all 
subsequent  pre-first-failure  increments  required.  All  post-first-failure 
applied  stress  increments  were  specified  as  500  lb/in.  . 

Total  failure  was  achieved  in  this  example;  the  resulting  composite 
stress-strain  curve  is  shown  in  Fig.  13.  As  in  computer  run  8A-10-17 
(see  Fig.  7),  the  plastic  component  of  the  computed  composite  strain  was 
in  error  by  a  factor  of  approximately  two.  This  accounts  for  the  small 
strains  in  Fig.  13,  since  the  curve  was  plotted  exactly  as  obtained  from 
the  computer  printout. 

As  in  8A-10-17 ,  first  failure  occurred  at  the  interface,  again  in 
element  139  (see  Fig.  3).  However,  the  applied  stress  at  first  failure, 
13,285  lb/in.2,  was  considerably  less  than  the  19,512  lb/in.  previously 
achieved,  as  was  the  composite  strain,  0.123  percent  versus  0.221  per¬ 
cent.  First  failure  occurred  during  applied  stress  increment  no.  9. 

Unlike  computer  run  8A-10-17,  the  crack  that  initiated  at  the  inter¬ 
face  immediately  propagated  along  the  interface,  upward  slightly  at  first, 
as  illustrated  in  Fig.  14.  However,  when  the  third  element  failed  during 
applied  stress  increment  no.  11,  a  total  of  nine  adjustment  increments 
were  required  before  equilibrium  was  again  achieved.  During  this  time 
the  crack  continued  to  propagate  along  the  interface,  primarily  down¬ 
ward;  the  condition  at  adjustment  increment  no.  11-6  is  shown  in  Fig.  15. 
Equilibrium  was  achieved  during  adjustment  increment  no.  11-9,  when  the 
crack  reached  the  lower  boundary,  as  indicated  in  Fig.  16.  However, 
during  the  next  two  applied  stress  increments,  the  crack  propagated  up¬ 
ward  along  the  interface  an  additional  two  elements,  causing  the  applied 
stress  to  drop  to  4563  lb/in.2.  This  marked  the  end  of  the  nearly  ver¬ 
tical  first  major  drop  in  applied  stress  beyond  first  failure. 

The  applied  stress  then  increased  during  the  next  three  applied 
stress  increments,  with  no  additional  elements  failing.  At  this  point 
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Fig.  13 —  Predicted  and  experimental  stress-strain  curves 
Computer  run  8A- 19-25 
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Fig.  14 — Failed  elements  and  octahedral  shear  stress  contours  (normalized) 
Computer  run  8A-19-25,  applied  stress  increment  11 
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Fig.  15 — Failed  elements  and  octahedral  shear  stress  contours  (normalized) 
Computer  run  8A- 19-25,  applied  stress  increment  11-6 
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Fig.  ]£,—  Failed  elements  and  octahedral  shear  stress  contours  (normalized) 
Computer  run  8A-19-25,  applied  stress  increment  11-9 
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the  crack  began  to  propagate  again,  moving  away  from  the  interface,  as 

shown  in  Fig.  17,  causing  another  drop  in  the  applied  stress.  During 

the  next  three  increments  the  applied  stress  again  increased,  with  no 

elements  failing.  But  the  stiffness  of  the  composite  was  considerably 

reduced,  as  shown  in  Fig.  13  by  the  slope  of  the  curve  segment  between 
2  2 

3043  lb/in.  and  3883  lb/in.  applied  stress.  This  was  due  to  the 
drastic  reduction  of  the  load-carrying  cross-sectional  area  caused  by 
the  extensive  crack  in  the  composite.  As  may  be  seen  in  Fig.  17,  only 
a  narrow  strip  of  material  between  the  crack  tip  and  the  upper  boundary 
remained  to  carry  the  applied  stress. 

During  the  final  three  applied  stress  increments,  nos.  24  through 
26,  the  crack  propagated  to  the  upper  boundary,  as  shown  in  Fig.  18. 

The  applied  stress  dropped  to  zero,  as  illustrated  in  Fig.  13,  defining 
total  failure.  As  Fig.  18  shows,  a  very  well-defined  crack  was  obtained. 

A  comparison  of  computer  runs  8A-10-17  and  8A-19-25,  which,  as 
previously  stated,  were  essentially  identical  except  for  the  magnitude 
of  the  applied  stress  increments  assumed,  reveals  several  points.  One 
is  that  a  basic  assumption  of  the  incremental  loading  procedure  is  that 
the  material  behavior  is  linear  within  each  applied  stress  increment. 

That  is,  the  slope  of  the  octahedral  shear  stress-plastic  octahedral 
shear  strain  curve  is  assumed  to  remain  constant  during  each  increment. 
This  effectively  means  that  each  material  element  follows  a  stress-strain 
curve  made  up  of  a  series  of  straight-line  segments  rather  than  the 
actual  input  curve.  A  new  straight-line  segment  of  proper  slope  is 
established  at  the  beginning  of  each  successive  applied  stress  increment. 
Thus,  large  increments  may  cause  the  element  material  response  to  oscil¬ 
late  excessively  about  the  smooth  input  stress-strain  curve.  During  a 
given  increment,  the  material  response  of  one  element  may  be  overshoot¬ 
ing  the  input  curve  while  an  adjacent  element  is  undershooting.  This 
results  in  a  nonuniform  local  stress  distribution,  particularly  in 
regions  where  the  stresses  are  on  a  portion  of  the  stress-strain  curve 
where  the  slope  is  rapidly  changing.  This  in  fact  did  occur  in  computer 
run  8A-19-25  as  first  failure  was  being  approached.  When  element  139 
failed  (at  a  lower  applied  stress  level  than  in  computer  run  8A-10-17 
because  of  the  increased  stress  due  to  an  overshoot)  and  the  required 
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Pigt  17 — Failed  elements  and  octahedral  shear  stress  contours  (  normalized) 
Computer  run  8A-  19-25,  applied  stress  increment  20-2 
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Fig.  16—  Failed  elements  and  octahedral  shear  stress  contours  ,  normal ized ) 
Computer  run  8A~  19-25,  applied  stress  increment  26-2  total  failure) 
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stress  adjustment  increment  was  performed,  the  stress  in  element  150 
overshot  the  actual  stress-strain  curve  excessively,  exceeding  the 
ultimate  stress  level.  This  element  failure  in  turn  caused  adjacent 
elements  to  fail,  setting  up  an  almost  self-propagating  effect.  This 
accounted  for  the  very  sudden  drop  from  the  maximum  applied  stress  down 
to  zero  stress  that  is  exhibited  in  Fig.  13. 

It  should  be  noted  that  a  similar  sharp  drop  occurred  in  Fig.  7 
when  the  crack  began  to  propagate  along  the  interface.  However,  in  that 
case,  considerable  strain  at  the  high  stress  levels  had  already  occurred 
while  the  crack  was  initially  propagating  in  the  matrix  region  away  from 


the  interface. 
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V.  DISCUSSION 


It  became  obvious  during  the  investigation  described  here  that 
much  additional  work  remains  to  be  done  in  problem  areas  that  range 
from  improving  the  accuracy  of  the  basic  finite  element  solution  tech¬ 
nique  to  developing  a  more  representative  model  of  the  propagating 
crack.  Many  of  these  areas  for  additional  research  are  well-defined 
and  require  only  the  time  to  develop  them.  The  solutions  to  other 
problems  are  not  so  clear. 

The  author  has  briefly  explored  some  of  these  research  areas, 
often  only  enough  to  discover  certain  inherent  problems.  Others  have 
been  outlined  in  detail  and  only  await  accomplishment.  The  following 
paragraphs  define  the  more  important  of  these  problems,  indicate  what 
has  been  done  to  date,  and  discuss  what  remains  to  be  done. 

FINITE  ELEMENT  METHODOLOGY 

As  discussed  in  detail  in  Sec.  II,  constant  strain  triangular  ele¬ 
ments  were  used  in  the  investigation.  Inaccuracies  in  the  incremental 
solution  results,  particularly  as  manifested  by  the  shear  strains, 
apparently  are  associated  with  the  finite  element  representation.  Higher 
order  elements  have  been  developed  that  more  closely  model  the  actual 
strain  variation  across  an  element.  The  linear  strain  triangular  ele¬ 
ment  has,  for  example,  received  considerable  attention.  Each  element 
involves  six  node  points,  one  at  the  midpoint  of  each  side  of  the 
triangle  and  one  at  each  apex.  The  author  has  not  attempted  to  use 
higher  order  elements  as  a  means  of  improving  accuracy,  and  this  should 
be  explored.  It  must  be  established  whether  using  additional  node  points 
in  a  higher  order  approximation  is  a  better  technique  than  simply  increas¬ 
ing  the  number  of  constant  strain  elements  proportionally.  A  thorough 

discussion  of  the  general  use  of  higher  order  elements  may  be  found  in 

(14) 

a  report  by  Felippa. 

Another  source  of  error  is  in  the  construction  of  the  grid  array 
itself.  Since  the  triangular  elements  can  be  assembled  in  an  arbitrary 
manner,  any  number  of  elements  can  have  a  common  node  point.  Thus, 
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various  patterns  are  created.  Walz ,  et  al .  have  investigated  this 

source  of  discretization  error  for  specific  grid  patterns.  A  brief 
study  was  also  performed  as  part  of  the  author’s  investigation.  Since 
inaccuracies  were  noted  in  even  the  first  (elastic)  applied  stress  in¬ 
crement,  only  elastic  behavior  was  analyzed. 

The  finite  element  equations  that  express  force  equilibrium  at  a 
given  node  point  in  terms  of  node  point  displacements  are  expanded  in 
Taylor  series  about  that  node  point  to  obtain  differential  equations. 
These  differential  equations  are  then  compared  with  the  classical  con¬ 
tinuum  equations  of  equilibrium.  Any  differences  represent  the  error 
terms . 

For  the  present  two-dimensional  plane  strain  analysis  and  an  arbi¬ 
trary  grid  pattern,  the  following  results  were  obtained: 


o 


o 


o 


o 


For  interior  node  points  surrounded  by  elements  all 


having  the  same  material  properties,  no  error  terms 
involving  first-order  derivatives  exist. 

For  interior  node  points  surrounded  by  elements  not 


all  having  the  same  material  properties,  e.g.,  at  a 


filament-matrix  interface,  there  are  error  terms  in 
volving  |^,  |^,  |p  and  |j,  which  contribute  to  errors 
in  ex,  ey,  and  Yxy. 

For  node  points  on  the  boundaries  x  =  0  and  x  -  a 

,  .  9u  ,  9v 

(see  Fig.  2)  ,  error  terms  involving  -g^*  and  -g—  exist 


in  the  equation  representing  equilibrium  in  the 

,  9u  ,  3v  . 

x-direction;  error  terms  involving  and  -g^  exist 
in  the  equation  representing  equilibrium  in  the 


y-direction. 

For  node  points  on  the  boundaries  y  -  0  and  y  =  b , 
there  are  error  terms  involving  -g~  and  -g^  in  the 

equation  representing  equilibrium  in  the  x-direction; 

9  9  v 

error  terms  involving  ^  and  exist  in  the  equation 
representing  equilibrium  in  the  y-direction. 


-77- 


Note  that,  except  for  certain  special  grid  patterns,  there  will  also 
be  error  terms  involving  second-order  derivatives.  All  of  these  error 
terms  are  independent  of  the  element  size — they  do  not  reduce  to  zero 
as  the  element  size  vanishes.  There  are  other  discretization  errors 
that  do  vanish  as  the  element  size  vanishes,  but  they  are  not  of  prin¬ 
cipal  concern  here. 

The  general  conclusion  is  that  additional  errors  are  introduced 
by  the  presence  of  material  interfaces  and  boundaries  (neither  of  which 
was  considered  in  Ref.  15)  and  cannot  be  avoided.  Their  severity  is 
dependent  upon  the  specific  local  grid  pattern,  but  they  are  difficult 
to  evaluate  analytically.  A  more  detailed  and  thorough  investigation 
of  this  aspect  of  the  accuracy  of  finite  element  approximations  and  the 
various  sources  of  error  is  needed. 

The  entire  problem  of  inherent  local  errors  leading  to  numerical 
instabilities  that  has  been  observed  in  this  study  is  by  no  means  unique. 
The  same  problem  has  been  reported  and  discussed  by  others,  including 
Zienkiewicz and  Argyris,^^  two  recognized  authorities  on  the  use 
of  finite  element  methods.  These  authors  were  concerned  primarily  with 
elastic  analyses,  for  which  only  one  matrix  inversion  is  needed.  The 
error  problem  is  compounded  when  an  incremental  solution  technique  in¬ 
volving  a  large  number  of  matrix  inversions  is  required,  as  in  the  appli¬ 
cation  of  the  tangent  modulus  method  of  elastoplastic  analysis  described 
in  this  report  (or  when  a  large  number  of  iteration  processes  are  required, 
as  in  the  method  of  initial  strains). 

Both  Zienkiewicz anc[  Argyris^^  have  suggested  the  use  of  some 
type  of  averaging  procedure  to  smooth  out  local  irregularities,  and  have 
indicated  successful  methods  for  doing  so.  Based  upon  the  general  sug¬ 
gestions  of  Argyris,  the  problem  might  be  approached  as  follows.  It  was 
noted  in  this  report  that  the  dominant  irregularities  occurred  in  the 
shear  strains,  particularly  in  regions  of  high  shear  strain  gradients. 

This  suggests  that  as  a  first  attempt,  only  the  shear  strains  for  each 
individual  boundary-value  problem  solution  be  averaged,  before  these 
solutions  are  combined  during  each  increment.  A  possible  procedure  is 


to: 
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1.  Assign  the  average  of  the  shear  strain  values  in  pairs  of 
adjacent  elements  to  the  midpoint  of  their  common  side.  At 
the  boundaries  of  the  material  region  being  analyzed,  this 
midpoint  value  will  equal  the  element  value,  because  of  the 
assumed  symmetry  across  the  boundary.  At  the  filament-matrix 
interface,  assign  the  element  value  to  the  midpoint  of  the 
side  along  the  interface,  i.e.,  do  not  average  across  the 
interface  since  the  materials  are  dissimilar  and  the  stresses 
discontinuous . 

2.  Average  the  values  on  the  three  sides  of  each  element  together 
and  assign  this  average  value  to  the  element. 

This  averaging  procedure  is  straightforward  and  can  be  readily  incorpor¬ 
ated  into  the  existing  computer  program.  Although  not  as  satisfying 
analytically  as  using  a  higher  order  element  representation,  it  may 
provide  numerical  results  of  satisfactory  accuracy  for  most  engineering 
applications  and  be  much  more  economical  of  computer  time  required. 

In  addition,  it  may  be  that  similar  irregularities  and  instabilities 
will  occur  if  higher  order  elements  are  introduced. 

CRACK  REPRESENTATIONS  AND  FAILURE  CRITERIA 

The  method  of  creation  of  a  crack  by  "failing"  an  element  and 
allowing  the  surrounding  material  to  adjust  to  the  failure  was  de¬ 
scribed  in  Sec.  II.  A  possible  alternative  method  of  modeling  a  crack 
is  discussed  in  Appendix  A.  In  some  actual  materials,  very  high  stress 
and  strain  gradients  exist  at  the  tip  of  a  sharp  crack.  The  adequacy 
of  using  finite  elements  to  model  these  gradients  should  be  explored 
further.  Obviously  the  "failed  element"  approach  used  in  the  study 
reported  here  does  not  represent  the  actual  crack  tip  geometry  very 
closely.  However,  in  most  real  materials  there  is  a  plastic  zone 
around  the  crack  tip  that  tends  to  reduce  the  magnitude  of  the  local 
stress  gradients,  and  may  make  it  an  adequate  model.  Certainly,  many 
practical  difficulties  can  arise  when  developing  a  crack  model ,  as 
pointed  out  in  Appendix  A. 
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The  definition  of  what  constitutes  material  failure  should  also  be 
explored  further.  Here,  exceeding  specified  ultimate  values  of  octa¬ 
hedral  shear  stress  and  octahedral  plastic  shear  strain  defines  failure. 
Other  theories  of  failure  have  been  proposed  and  they  should  be  evalu¬ 
ated  for  the  application  discussed  in  this  report.  It  is  hoped  that 
the  sharp  increase  in  interest  in  fracture  mechanics  concepts  also  will 
prompt  applications  relevant  to  this  problem. 


GENERALIZED  PLANE  STRAIN 


This  analysis  assumes  a  condition  of  plane  strain,  i.e.,  that  the 

normal  strain  e  in  the  z-direction  (the  direction  of  the  filament  re- 
z 

inforcement)  is  zero  throughout  the  composite  body.  Correspondingly, 
the  z-component  of  displacement,  w,  is  zero  everywhere.  When  transverse 
normal  loadings  and  a are  applied,  however,  O z  components  of  stress 
are  induced.  Thus,  an  average  stress  in  the  z-direction,  a  ,  exists. 

In  actual  applications  it  may  be  desirable  to  keep  this  average 

stress  equal  to  zero  to  simulate  a  real  loading  condition.  This 

can  be  done  by  assuming  a  uniform  nonzero  axial  displacement  w  rather 

than  a  zero  value.  Such  a  condition  is  referred  to  as  generalized 

plane  strain.  The  magnitude  of  w  required  to  achieve  O  =  0  is  ini- 

tially  unknown  and  must  be  solved  for.  For  the  problem  considered 

here,  this  will  require  the  solution  of  an  additional  boundary-value 

problem  for  each  increment,  with  Aw  an  arbitrarily  specified  constant 

and  the  normal  displacement  increments  Au  and  Av  along  the  boundaries 

x  =  a  and  y  =  b,  respectively,  set  equal  to  zero.  The  constitutive 

relations  must  also  be  modified  slightly  to  accommodate  the  nonzero 

value  of  £  .  A  more  detailed  discussion  is  contained  in  Ref.  11. 
z 

The  additional  boundary-value  problem  is  then  combined  with  the 
appropriate  two  problems  defined  in  Sec.  II  (Problems  1  and  2  for  an 
applied  stress  increment,  or  Problems  2  and  3  for  an  adjustment  incre¬ 


ment)  to  obtain  the  desired  values  of  Aa^  and  Aa^ ,  and  the  generalized 
plane  strain  condition,  A =  0. 

Generalized  plane  strain  was  not  used  in  the  present  analysis  be¬ 


cause  the  required  solution  of  three  rather  than  two  boundary-value 


problems  per  applied  stress  or  adjustment  increment  would  have  increased 
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the  computation  time  by  approximately  fifty  percent.  Interest  here 
has  been  concentrated  on  developing  the  concept  of  crack  propagation 
rather  than  on  producing  an  all-encompassing  solution  technique.  Gen¬ 
eralized  plane  strain  will  be  a  logical  extension,  however,  particularly 
because  of  the  interest  in  combined  loadings,  as  discussed  below. 


COMBINED  LOADINGS 

In  a  laminated  composite  system,  the  individual  laminae  are  usually 
subjected  to  a  combination  of  axial,  transverse,  and  shear  loadings  due 
to  interactions  between  the  laminae,  even  when  the  composite  body  is 
subjected  to  a  simple  loading  state. 

The  same  basic  methodology  developed  here  for  crack  propagation  in 
a  composite  subjected  to  transverse  normal  loading  can  be  applied  to 
combined  loading.  The  generalized  plane  strain  condition  discussed  in 
the  previous  subsection  is,  in  fact,  one  special  case  of  axial  loading, 
i.e.,  when  0^  =  0 .  Thus,  by  specifying  =  0^  =  0  and  0^  4-  0 ,  axial 
loading  of  a  unidirectional  composite  is  obtained.  By  specifying  arbi¬ 
trary  values  of  0  ,  0  ,  and  0  ,  combined  axial  and  transverse  normal 

x  y  z 

loading  is  obtained.  The  longitudinal  shear  loadings  T  and  T  re- 

z  x  z  y 

main  to  be  added. 

A  finite  element  analysis  of  the  elastoplastic  behavior  of  a  uni¬ 
directional  composite  subjected  to  longitudinal  shear  loading  was  de¬ 
scribed  in  Ref.  11  using  the  method  of  initial  strains,  and  in  unpublished 

k 

work  by  Adams  and  Wilson  using  the  tangent  modulus  method.  Combined 
axial,  transverse,  and  shear  loadings  of  an  elastoplastic  material  were 
also  considered  in  Ref.  11. 


These  same  methodologies  can  be  applied  to  the  problem  of  crack 
initiation  and  propagation.  The  addition  of  specified  increments  of 
applied  longitudinal  shear  stress  components  to  the  combined  loading  will 
require  the  solution  of  one  additional  boundary-value  problem  per  solu¬ 


tion  increment  if  only  one  of  the  two  components  At ■  or 
included,  and  two  additional  solutions  if  both  are  to  be 


Ax  is  to  be 
zy 

included.  Thus, 


k 

D.  F.  Adams  and  H.  B.  Wilson,  Jr.,  "Inelastic  Analysis  of  a  Uni¬ 
directional  Composite-Longitudinal  Shear  Loading,"  unpublished  work, 
Aeronutronic  Division,  Philco-Ford  Corporation,  August  1966. 
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a  complete  combined  loading  problem  involving  the  applied  stress  incre¬ 
ments  A o9  Aa  ,  Aa  ,  Ax  ,  and  Ax  will  require  the  solution  and  com- 
x  y  z  /jX  zy 

bination  of  five  separate  boundary-value  problems  per  solution  increment. 
However,  longitudinal  shear  loading  involves  only  one  unknown  displace¬ 
ment  component  per  node  point  (a  displacement  in  the  z-direction)  rather 
than  two  as  in  the  case  of  axial  or  transverse  normal  loading  (displace¬ 
ments  in  both  the  x-  and  y-directions) .  Thus,  only  an  N  x  N  matrix 
rather  than  a  2N  x  2N  matrix  needs  to  be  inverted  to  obtain  a  solution 
(N  represents  the  total  number  of  node  points  in  the  finite  element 
array) .  This  suggests  that  the  computation  time  required  to  solve  a 
longitudinal  shear  boundary-value  problem  will  be  only  about  one-fourth 
as  long , 

The  general  expressions  for  octahedral  shear  stress  and  octahedral 
plastic  shear  strain  given  by  Eqs .  (7)  and  (6),  respectively,  are  used 
to  define  material  yield  and  failure.  The  consideration  of  crack  ini¬ 
tiation  and  propagation  under  combined  loading  should  be  a  straight¬ 
forward  extension  of  the  investigation  given  here  and  that  found  in 
Ref.  11. 

THERMAL  EFFECTS 

Thermal  effects  are  of  special  interest  in  studies  of  composite 
material  micromechanical  behavior  because  the  differences  in  coefficients 
of  thermal  expansion  among  the  constituents  almost  always  result  in  local 
stresses  being  introduced  during  a  temperature  change,  even  a  uniform 
temperature  change.  Most  composite  materials,  both  metal  and  polymer 
matrix  composites,  are  produced  at  elevated  temperatures  and  subsequently 
are  used  in  room  temperature  applications.  Thus,  even  the  as-received 
composite  material  contains  (thermal)  residual  stresses.  If  the  com¬ 
posite  is  also  subjected  to  thermal  environments  during  service,  addi¬ 
tional  thermal  stresses  are  developed. 

Since  the  thermal  stresses  introduced  in  fabrication  are  developed 
while  the  composite  is  not  being  subjected  to  significant  mechanical 
loads,  i.e.,  during  cool-down  from  fabrication  temperatures,  the  con¬ 
stituent  materials  will,  at  least  initially,  exhibit  linear  elastic 
behavior.  In  most  instances,  these  thermal  stresses  will  not  be  high 
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enough  to  cause  yielding.  This  is  particularly  true  because  the  elastic 
stiffness  (Young !s  modulus)  of  the  constituents,  especially  the  matrix 
material  that  must  flow  during  the  compositing  process,  will  be  lower 
at  the  fabrication  temperature.  Thus,  a  linearly  elastic  thermostress 
analysis  is  often  applicable.  The  variation  of  material  properties  with 
temperature,  particularly  the  Young ?s  modulus,  should  be  taken  into 
account,  however,  which  complicates  the  analysis. 

A  linearly  elastic  micromechanical  analysis  of  thermal  residual 
stresses  was  included  in  Ref.  2.  The  variation  of  the  Young  s  modulus 
with  termperature  was  taken  into  account  in  an  approximate  way  by  using 
room  temperature  material  properties  but  an  assumed  temperature  change 
less  than  the  actual  value.  The  actual  temperature  change  associated 
with  cool-down  can  be  applied  in  increments  using  the  analysis  described 
here,  with  the  material  properties  modified  at  the  beginning  of  each 
increment.  This  procedure  will  provide  a  more  accurate  representation 
of  the  actual  physical  behavior. 

Similarly,  thermal  loading  need  not  be  limited  to  elastic  response. 
By  suitably  modifying  the  elastoplastic  analysis,  arbitrary  (including 
nonuniform)  temperature  variations  can  be  combined  with  the  mechanical 
loading  increments.  This  will  require  the  specification  of  appropriate 
material  properties,  including  a  complete  stress-strain  curve,  for  each 
temperature  increment  used.  Since,  in  the  tangent  modulus  method,  the 
stiffness  matrix  is  reconstructed  and  inverted  in  each  increment  anyway, 
this  modification  will  not  involve  excessive  additional  computational 
time.  However,  the  inclusion  of  a  thermal  loading  (temperature  change) 
during  a  solution  increment,  just  as  any  of  the  other  loadings  pre¬ 
viously  discussed,  will  require  the  solution  of  an  additional  boundary- 
value  problem,  to  be  combined  with  the  others.  The  conditions  for  the 

problem  will  be:  zero  shear  stress  T  and  zero  normal  displacement 

xy 

along  all  four  rectangular  boundaries,  and  specification  of  the  desired 
incremental  temperature  change  in  each  element. 

The  governing  field  equations  must  also  be  modified  for  this 
boundary— value  problem  to  include  the  thermal  strain  effect.  For  ex¬ 
ample,  the  strain-displacement  relation  given  by  Eq .  (35)  must  be 
amended  to  include  a  thermal  strain  term  [a]AT,  i.e., 
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IeJ  -  tOJ  [63  +  MAT 

where 

[a]  - 

represents  the  coefficient  of  thermal  expansion  and  AT  is  the  incremental 
temperature  change.  The  constitutive  relation,  Eq .  (30),  remains  un¬ 
changed.  Following  through  the  analysis  of  Sec.  II,  Eq .  (42)  becomes 

[F]  =  [K]  [6]  +  An[6]T  [H]  [a]  AT 

The  additional  term  is  a  known  constant  for  each  element  during  the 
increment.  Thus,  no  special  complications  will  occur  in  formulating  a 
solution  procedure. 

DIAMOND  PACKING  ARRAYS 

Filament  packing  geometries  and  the  desirability  of  maintaining  a 
regular  periodic  array  was  discussed  in  Sec.  II.  Although  this  investi¬ 
gation  has  been  limited  to  rectangular  arrays,  diamond  arrays,  of  which 
a  hexagonal  array  is  a  special  case,  are  of  equal  interest.  As  noted 
on  the  left  side  of  Fig.  1c  and  in  the  text  of  Sec.  II,  this  analysis 
also  can  be  used  to  treat  diamond  arrays.  It  is  possible  to  reduce  the 
area  in  the  first  quadrant  that  must  be  analyzed  to  one-half  by  redefining 
the  boundary-value  problem  to  include  only  that  region  indicated  by  the 
solid  lines  on  the  right  side  of  Fig.  lc.  This  region  of  the  first 
quadrant  is  redrawn  in  Fig.  19. 

The  greater  difficulty  in  analyzing  a  diamond  array  as  compared 
with  a  rectangular  array  is  associated  with  the  fact  that  the  boundary 
represented  by  the  line  COD  in  Fig.  19  does  not  remain  a  straight  line 
under  loading.  However,  it  does  deform  antisymmetrically  about  the  mid¬ 
point  0,  i.e.,  the  displacements  of  two  points  on  line' COD  at  equal  dis¬ 
tances  on  opposite  sides  of  point  0  will  be  equal  in  magnitude  but 
opposite  in  direction.  Also,  because  of  the  assumed  geometric  symmetry, 


a 

a 

o 
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forces  acting  on  these  points  will  be  of  equal  magnitude  and  direction. 
Thus,  by  constructing  node  points  on  boundary  COD  symmetric  about  the 
midpoint  0,  both  displacement  and  force  relations  between  pairs  of  node 
points  can  be  established. 

This  is  adequate  information  to  solve  the  boundary-value  problem. 
However,  there  are  also  certain  practical  difficulties.  These  relation¬ 
ships  are  of  a  different  form  than  those  making  up  the  remainder  of  the 
stiffness  matrix,  and  thus  the  solution  procedure  must  be  modified.  Node 
points  as  widely  separated  as  C  and  D  are  now  related,  a  fact  to  be  con¬ 
sidered  when  numbering  the  node  points  so  as  to  minimize  the  matrix 
bandwidth . 

An  alternative  scheme,  which  eliminates  the  need  for  including  dis¬ 
placement  and  force  equivalence  relations  in  the  stiffness  matrix  di¬ 
rectly,  is  to  add  additional  node  points  and  elements  on  the  outside  of 
boundary  COD.  The  node  points  on  COD  are  still  placed  symmetrically 
about  midpoint  0.  Additional  fictitious  node  points,  symmetric  about 
point  0  to  node  points  inside  the  material  region,  are  added  outside  so 
that  elements  can  be  constructed  that  completely  surround  the  node  points 
on  one-half  of  boundary  COD,  say  between  C  and  0,  as  shown  in  Fig.  19. 

The  displacements  of  the  boundary  node  points  between  0  and  D  are  re¬ 
lated  to  those  of  their  symmetric  counterparts  between  C  and  0  as  be¬ 
fore — they  are  of  equal  magnitude  but  opposite  direction  relative  to 
point  0.  Equilibrium  equations  can  then  be  written  for  the  node  points 
on  the  boundary  between  C  and  0  just  as  if  they  were  interior  nodes,  and 
the  displacements  of  the  fictitious  nodes  are  eliminated  by  using  the 
known  relationships  to  their  symmetric  counterparts  in  the  interior. 

This  scheme  eliminates  the  need  for  satisfying  force  equivalences 
across  the  midpoint  of  the  boundary  COD;  the  resulting  stiffness  matrix 
retains  a  form  similar  to  that  for  a  rectangular  array.  Numbering  the 
node  points  to  keep  the  matrix  bandwidth  small  remains  an  important  con¬ 
sideration,  just  as  it  is  for  the  first  scheme. 

OTHER  TOPICS 


There  are  a  number  of  other  topics  requiring  additional  investiga¬ 
tion,  only  a  representative  sample  of  which  will  be  mentioned  here.  For 


-86- 


example,  this  investigation  has  only  considered  regular  arrays  of  fila¬ 
ments,  whereas  in  actual  composite  materials  the  filaments  are  often 
not  regularly  spaced.  The  influence  of  random  filament  packing  on 
elastic  stiffness  properties  was  investigated  in  Ref.  3.  The  influence 
of  randomness  on  crack  propagation  and  composite  ultimate  strength 
should  also  be  studied. 

Anisotropy  of  the  constituents  is  also  a  consideration.  Graphite 
filaments,  for  example,  are  known  to  be  anisotropic  due  to  their  crystal¬ 
line  structure. 

Certain  filament  reinforcements,  e.g.,  graphite  produced  from  a 
rayon  precurser,  have  noncircular  and  often  very  irregular  cross-sectional 
shapes.  The  analysis  described  in  this  report  can  readily  handle  an  ir¬ 
regular  geometry.  A  study  should  be  performed  to  determine  the  influence 
of  an  irregular  shape  on  crack  propagation  and  ultimate  strength. 

Voids  in  the  matrix  -material  and  local  unbonded  regions  at  the  fila¬ 
ment-matrix  interface  are  examples  of  defects  often  present  in  an  actual 
composite.  A  systematic  study  of  their  influence  on  crack  propagation 
and  ultimate  strength  should  be  conducted. 
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Appendix  A 

INVERSE  OF  THE  CONSTITUTIVE  EQUATION 


The  inverse  of  Eq .  (25)  can  be  obtained  by  first  expanding  it  into 

•  •  i 

the  three  equations  it  represents,  G22*  anc*  el2*  anc*  t^en  s°lving 

these  equations  for  0^,  ®  22s  anc*  °12  successive  elimination.  Be¬ 
cause  such  a  process  is  lengthy,  it  is  much  simpler  to  invert  exclu¬ 
sively  in  index  notation,  as  follows.  Repeating  Eq .  (25), 


1+V 


'aB 


(a 


a8 


-  wyy5o6) 


1 

+  —  t  at  -a  . 

D  a8  yS  yS 


from  which  the  following  two  equations  are  obtained: 


aa 


=  C 


1  +  v 


)(l  -  2v)a  + 


aa 


t  t  r 

aa  y  6 


and 


(A-l) 


•  l~f*V  *  •  i  • 

fca8ea8  -  E  ^ta6aa8  vtaaaYY^  +  D  ta8ta8tY<ScrY5 


(A-2) 


Letting  x  =  and  y  =  »  Eqs .  (25),  (A-l),  and  (A-2)  can  be 

written  as 


SB  "  TT  (°aB  ‘  ^aB^  +  S  taBy 


(A-3) 


aa 


(1  +  v)(l  -  2v)  .  1  , 

- - -  x  +  —  t  y 

E  D  aaJ 


(A-4) 


t  a£  0  =  (y  -  vt  x)  +  7T  t  at  ay 

a8  a8  E  aa  D  a8  a8 


(A-5) 


Eliminating  x  and  y,  solving  for  and  simplifying  give  the  desired 

inverse,  Eq.  (26). 
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First 


Alternatively 
,  determine  s 


,  the  inverse  of  Eq . 

cB'aB  £ro"  Eq'  <25)’ 


(25)  can  be  obtained  as  follows. 


sa6eotB 


1-t-V 

E 


.saB°oiB 


-  vs  no  6  n 
a8  YY  a6. 


+  D  ta6tY($SoiB0Y(5 


noting  that,  as  usual,  the  repeated  indices  indicate  summation.  Thus, 
the  second  term  in  the  brackets  can  be  expressed  as  v(s  11  +  S22)CYY’  °r 

as  ~Vs33°a86aB’  since  (sll  +  S22}  =  ~s33  and  °YY  7  °a86a6-  The  exPres" 
sion  within  the  brackets  becomes  [  (s^  +  Vs 33^0,3^ aot6 ^  '  But  the  cJuantity 

within  the  parentheses  here  is  equal  to  t^g ,  as  previously  defined. 

Therefore,  the  full  equation  can  be  rewritten  as 


saBeaB 


l+v 


taBSaB ' 

D 


taBaaB 


Setting  the  terms  in  parentheses  over  the  common  denominat  or  l+v  D  and 
expanding  t^s^  gives 


saBea3 


A  +  Es 


2 

33 


+ 


E 

1+V 


+  2s 


2 

12 


+  VS 


33 


taB°otB 


l+v 


o 

where  the  term  s^0  has  been  added  and  subtracted  within  the  brackets. 

33  2 

The  first  four  terms  within  the  brackets  are  now  equal  to  3Tq.  Also, 

since  (s,-,  +  s00)  -  -sQQ,  the  remaining  terms  within  the  brackets  re- 
11  22  ~  33  £ 

duce  to  -(1  +  v)s^.  This  term,  when  multiplied  by  the  yyy  appearing 
outside  the  brackets,  becomes  -Es^,  which  cancels  with  the  second 
term  within  the  braces.  Thus,  the  complete  equation  reduces  to 


s  D  -  (A  + 

aB  aB 


l  +  v 


-2*  taBaaB 
l+v 


which  can  be  rewritten  as 
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ta8aaB _ sa8£a8 _ 

D  Q  2  ,  ,  .A 

3T0  +  (1  +  V)  E 


(A-6) 


This  equation  relates  the  stress  a  *  in  Eq.  (25)  to  the  strain  £  ^ , 
The  next  step  is  to  relate  a  in  Eq.  (25)  to  £yy  *  Writing 

Eq .  (25)  for  £  , 

H  aor 


e  =  ~  (a  -  va  6  )  +  t  „ 

aa  E  aa  YY  aa  D  aa 


Simplifying  the  first  term  and  substituting  Eq .  (A-6)  into  the  second 

gives 

*  1+v  , ,  _  N  *  .  sy6£Y|S _  . 

£aa  E  V  aaa  ,  2  n  .A  aa 

3T0  +  (1  +  V)  E 


Substituting  the  expansion  of  t^  below. 


taa  ^sll  +  Vs33*  +  (s22  +  Vs33^ 


(sll  +  s22}  +  2Vs33 


s33  +  2vs33 


■(1  -  2v)s, 


into  the  above  equation  and  solving  for  a '  gives 


=  E _  •  (1  2v)s33s^6£y<s 

aa  (1  +  v)(l  -  2v)  eaa  ^2  +  ^  +  ^  A 

•  0  E 


(A- 7 ) 


Substituting  Eqs .  (A-6)  and  (A-7)  into  Eq .  (25)  gives  . 
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•  1+v  ! 

(  •  VE 

(1  -  2v)s33sY6eY61 

6  „ 

ea8  E 

|aa8  “  (1  +  v)(l  -  2v) 

^  +  3Tg  +  (1  +  v)  |  _ 

aB 

+  t 


sy6£y6 


010  3Tq  +  (1  +  v)  | 


Solving  for  0 


a  8 


aB  1+v 


v 


(Vs336a6  '  taB)8Y6£Yg 


ea6  +  1-2V  Yy5o6  +  3t2+(i  +  v)|  j 


Since  the  quantity  in  parentheses  in  the  numerator  of  the  last  term  is, 
by  definition,  -s^g,  the  equation  can  be  written  as 


a 


aB  1+v 


e  n  + 


v 


£  <5 


sa8SY6£Y(S 


'aB  1  1-2V  YY  3  2  +  (1  +  v)  A 

0  ti 


(A-8) 


which  is  the  inverse  of  Eq.  (25). 
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Appendix  B 

ALTERNATIVE  METHODS  OF  MODELING  A  CRACK 


In  addition  to  the  "failed  element"  concept  of  modeling  a  local 
failure,  which  was  used  in  the  investigation  discussed  here,  other 
methods  were  considered,  one  of  which  will  be  briefly  outlined  below. 
This  discussion  is  included  to  point  out  some  of  the  practical  diffi¬ 
culties  involved  in  modeling  a  crack  and  to  serve  as  a  starting  point 
for  additional  investigations  by  the  interested  reader.  There  is 
clearly  a  need  for  much  more  study  in  this  problem  area. 

When  the  stresses  in  a  given  element  exceed  the  ultimate  strength 
of  the  material,  rather  than  dropping  out  that  element  to  create  a 
crack,  an  additional  node  point  can  be  added  at  each  of  two  of  the 
three  corners  of  the  triangular  element.  These  node  points  are  co¬ 
located  with,  but  not  attached  to,  the  existing  ones.  At  each  of  these 
two  locations,  one  node  point  is  then  associated  with  certain  of  the 
surrounding  elements,  and  the  other  node  point  with  the  remainder. 

This  is  illustrated  in  Fig.  B-l,  which  represents  a  local  region  of 
the  full  grid  shown  in  Fig.  3.  Suppose  that  failure  is  indicated  in 
element  163.  Additional  node  points  177  and  178  are  added  at  nodes  68 
and  78,  respectively.  (See  Fig.  B-l(b) ;  the  nodes  of  the  pairs  68,  177 
and  78,  178  are  shown  separated  from  each  other  for  clarity,  and  the 
shaded  area  represents  the  crack.)  In  this  representation,  a  crack 
of  zero  initial  width  is  created  and  no  material  is  removed. 

However,  there  are  certain  procedural  difficulties  created.  For 
indicated  failure  in  element  163,  one  of  the  other  two  sides  could  have 
been  selected  as  the  crack  site  instead,  resulting,  for  example,  in  a 
crack  through  nodes  86,  77,  68,  59  or  nodes  76,  77,  78,  79  rather  than 
nodes  88,  78,  68,  58.  These  differences  are  significant,  resulting  in 
possible  cracks  almost  perpendicular  to  the  one  indicated  in  Fig.  B-l(b). 

These  are  not  the  only  three  crack  path  choices  possible,  although 
they  are  perhaps  among  the  more  obvious  for  the  particular  local  region 
selected  as  an  example  here.  For  instance,  the  crack  between  nodes  78 
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and  68  could  have  been  assumed  to  run  from  node  68  to  node  67,  59,  or 
69  rather  than  to  node  58,  Similar  choices  were  possible  for  the  other 
end.  Thus,  a  scheme  would  have  to  be  developed  for  selecting  the  most 
logical  path,  perhaps  based  upon  the  stress  gradients  existing  around 
the  element  in  which  failure  is  indicated. 

Figure  B-l  demonstrates  another  potential  problem.  An  assumed 
crack  along  either  of  the  other  two  sides  of  element  163  would  neces¬ 
sarily  extend  either  into  the  adjacent  material  region  (from  matrix 
into  filament)  or  along  the  interface.  In  an  actual  material  this  may 
not  be  realistic,  particularly  if  one  material  has  significantly  dif¬ 
ferent  properties  than  the  adjacent  one,  as  is  usually  the  case.  The 
problem  stems  from  the  fact  that,  although  failure  is  indicated  in  only 
one  element,  the  crack  is  modeled  to  extend  along  the  sides  of  as  many 
as  three  adjacent  elements,  as  may  be  seen  in  Fig.  B-l(b).  To  use 
this  model,  a  method  would  also  have  to  be  devised  to  insure  that  the 
crack  created  by  colocated  node  points  can  only  open  up;  the  material 
region  on  one  side  of  the  crack  cannot  be  allowed  to  move  into  the 
material  region  on  the  other  side,  based  upon  physical  considerations. 

A  serious  computational  difficulty  also  occurs  if  the  "double  node" 
crack  representation  is  used.  As  was  discussed  in  Sec.  II,  the  band¬ 
width  of  the  diagonalized  stiffness  matrix  [K-^]  in  Eq .  (51)  must  be 
kept  to  a  minimum  to  achieve  an  efficient  solution.  This  bandwidth  is 
proportional  to  the  maximum  difference  between  the  highest  and  lowest 
node  point  number  associated  with  each  of  the  various  elements  in  the 
array.  For  the  grid  array  shown  in  Fig.  3,  this  maximum  difference  is 
only  11.  However,  if  node  points  177  and  178  were  added  as  indicated 
in  Fig.  B-l(b)  ,  this  maximum  difference  would  increase  to  119  (for 
element  151).  Thus,  some  scheme  would  have  to  be  developed  for  auto¬ 
matically  renumbering  all  of  the  node  points  each  time  additional  node 
points  were  added.  At  present,  not  even  the  initial  set  of  node  points 
are  numbered  automatically;  because  of  the  computational  problems  in¬ 
volved,  they  are  numbered  by  hand. 

An  alternative  scheme  would  be  to  double-number  (consecutively)  all 
of  the  node  points,  and  use,  for  example,  only  the  odd  numbers  initially. 
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The  even-numbered  nodes  could  be  introduced  as  needed  to  model  the 
crack  formation.  This  scheme  would  approximately  double  the  matrix 
bandwidth,  however,  which  would  make  it  of  questionable  utility. 
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